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1 .0  INTRODUCTION 


To  obtain  useful  Images  of  satellites  from  the  ground,  large 
apertures  are  needed.  According  to  the  laws  of  diffraction,  to  obtain 
a  resolution  p  at  a  range-to-target  R  with  light  of  wavelength  X,  an 
aperture  of  diameter  about  D  «  XR Ip  is  required,  At  optical  or  near- 
infrared  wavelengths,  we  would  therefore  ,  need  an  aperture  on  the  order 
of  one  meter  In  diameter  for  low-altitude  earth-orbiting  satellite,  and 
an  aperture  the  size  of  a  football  field  for  geosynchronous  satellites. 
Conventional  telescopes  fall  to  achieve  this  for  two  reasons.  First, 
for  the  case  of  high-altitude  satellites,  the  required  large  apertures 
are  well  beyond  the  current  state  of  the  art.  Second,  for  all  cases, 
atmospheric  turbulence  limits  the  resolution  to  an  effective  aperture 
diameter  of  rQ  (Fried's  parameter),  which  is  typically  in  the  range  of 
0.05  to  0.20  meters,  a  small  fraction  of  what  is  needed.  Compensated 
imaging  systems,  consisting  of  wavefront  sensors  coupled  with  adaptive 
optics  to  correct  for  atmospheric  turbulence  in  real-time,  may  work 
well  for  the  low-altitude  case,  but  that  technology  does  not  scale  well 
for  the  high-altitude  case. 

This  report  describes  new  methods  we  have  developed  for 
reconstructing  fine-resolution  images  of  satellites  which  circumvent 
both  the  problems  of  atmospheric  turbulence  and  obtaining  large 
apertures  with  today's  technology.  These  methods  are  also  advantageous 
because  they  employ  the  simplest,  least  expensive  receiver  possible. 
Hardware  complexity  is  minimised,  however,  at  the  expense  of  software 
complexity  and  computing  requirements.  This  trade-off  is  increasingly 
advantageous  as  computers  have,  been  evolving  much  faster  than  optics 
and  detectors.  We  have  also  developed  methods  that  take  data  from 
sensors  based  on  other  imaging  concepts  und  reduce  the  phase  errors 
that  are  present  in  the  data. 
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While,  In  this  report,  we  concentrate  on  the  problem  of  imaging 
satellites,  the  methods  described  here  can  also  be  used  for  other 

i  i  ■  ■ 

applications  such  as  SOI  discrimination  and  tactical  Imaging. 

The  assumed  optical  system  for  the  majority  of  this  report  Is  as 
follows.  A  pulsed  laser  of  long  coherence  length  Illuminates  the 
target  object.  The  Intensity  ’  of  the  backscattered  radiation  Is 
detected  by  an  array  of  simple  light-bucket  detectors  on  the  ground, 
distributed  over  an  area  of  diameter,  D.  Alternatively,  for  the  low- 
altitude  case,  the  radiation  can  be  collected  by  a  telescope,  but  the 
detector  array  Is  placed  In  a  plane  conjugate  to  (l.e. ,  at  an  Image  of) 

the  aperture  plane,  rather  than  In  the  usual  focal  plane. 

.  «  . 

Since  the  phase  errors  due  to  atmospheric  turbulence  are  Introduced 
In  a  volume  relatively  near  to  the  aperture  plane,  the  detected 
Intensities  are,  to  first  order,  unaffected  by  the  atmospheric 

turbulence.  (The  Intensity  Is  affected  only  to  the  extent  that 

■  •  ,  .  ■  .  .  'i 

scintillation,  caused  by  strong  upper-atmospheric  turbulence,  Is 
present.)  If  the  phase  associated  with  the  Intensities  could  be 
retrieved,  then  by  digitally  back-propagating  the  wavefront  at  the 
aperture  plane  to  a  plane  at  the  target  object  (essentially  a  Fourier 
transform),  we  could  obtain  a  diffraction-limited  Image.  This  Image 
would  be  complex-valued  and  would  suffer  from  speckle.  Since  the 
realization  of  the  Image  speckle  pattern  would  change  for  each  laser 
pulse  (as  the  object  rotates  slightly  or  translates),  by  averaging  over 
the  Intensities  of  several  such  Images  we  can  average  out  the  speckles 
and  obtain  the  equivalent  of  an  Incoherent,  speckle-free  Image  of  the 
object.  The  required  phases  can  be  computed  using  one  of  the  phase- 
retrieval  (Image  reconstruction)  algorithms  developed  under  this 
effort.  They  require  either  a  glint  or  glints  to  be  present  on  the 
object  or  to  have  measured  partial  Information  about  the  phase  of  the 
wavefront  In  the  aperture  plane. 
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In  a  second  Imaging  mode,  the  multiple  realizations  of  the 
aperture-plane  Intensities  can  be  averaged  so  as  to  obtain  the  modulus 
(magnitude)  of  the  Fourier  transform  of  the  Incoherent  Image.  Again, 
by  the  phase  retrieval  algorithm,  the  phase  of  the  Fourier  transform 
can  be  retrieved,  allowing  the  reconstruction  of  a  diffraction-limited 
Incoherent  Image.  We  have  called  this  latter  mode  "Imaging 
correlography."  It  has  the  advantage  of  working  under  much  broader 
circumstances  (no  glints  or  partial  phase  are  required),  but  has  the 
disadvantage  of  requiring  a  large  number  of  laser  pulses  to  build  up  a 
sufficient  slgnal-to-nolse  ratio. 

This  report  describes  the  techniques  developed  for  these  two  novel 
Imaging  modes  and  phase-error  correction  techniques  that  can  be  applied 
to  other  Imaging  sensors.  It  also  describes  several  associated  topics, 
Including  the  comparison  of  several  competing  Imaging  approaches. 

In  overview,  the  Imaging  approaches  we  have  developed,  which 
utilize  phase  retrieval  algorithms,  make  possible  the  reconstruction  of 
fine-resolution  Images  of  satellites  using  hardware  technology  that  Is 
available  today.  The  approaches  require  only  light-bucket  detectors 
which  require  no  phasing  and  do  not  have  to  be  very  fast.  The  cost  of 
such  a  system  would  be  a  small  fraction  of  the  cost  of  most  competing 
approaches.  Alternatively,  Improved  Images  can  be  obtained  from  other 
imaging  sensors  by  correction  of  residual  phase  errors. 

Section  2  of  this  report  contains  a  brief  summary  of  the  research 
accomplishments.  Individual  topics  are  described  In  detail  In  Sections 
3  to  8.  References  are  found  at  the  end  of  each  section. 
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2.0  SUMMARY  OF  ACCOMPLISHMENTS 


In  this  section,  we  briefly  summarize  the  accomplishments  of  our 
research  on  active  Imaging.  Details  are  given  In  the  sections  that 
follow. 

As  described  In  Section  ,3,  we  developed  a  new  Imaging  modality 
called  Imaging  correlography.  In  It,  multiple  arrays  of  aperture-plane 
Intensity  measurements  of  laser  pulses  backscattered  by  the  object  are 
collected.  By  averaging  over  the  autocovariances  of  these  Intensity 
measurements,  we  arrive  at  an  estimate  of  the  squared  modulus  of  the 
Fourier  transform  of  the  Incoherent  object  (the  object  reflectivity 
function  for  incoherent  Illumination).  From  these  data,  a  fine- 
resolution  nonspeckled  Image  can  be  reconstructed  using  the  Iterative 
Fourier  transform  (phase  retrieval)  algorithm. 

The  bas»1c  theory  of  Imaging  correlography  was  developed.  Averaging 
In  both  the  aperture  plane  and  the  Fourier  transform  of  the  aperture 
plane  were  analyzed.  Multiple  realizations  of  complex-valued,  laser- 
illuminated  reflectivity  functions  of  a  satellite  model  were  computer 
simulated,  and  the  aperture-plane  data  were  simulated.  The  averages 
were  computed,  and  Wiener  filtering  of  the  resultant  Fourier  modulus 
estimates  was  performed.  Images  were  reconstructed  from  these  data  for 
various  numbers  of  frames  (laser  pulses).  When  a  large  number  of 
frames  were  processed,  high-quality,  fine-resolution  Images  were 
successfully  reconstructed. 

The  same  experiments  were  repeated  for  the  case  of  a  sparse 
collecting  array  consisting  of  a  Golay  arrangement  of  six  subapertures. 
The  ability  of  Wiener  filtering  to  correct  for  the  effects  of  the  MTF 
of  the  sparse  array  and  the  ability  to  reconstruct  fine-resolution 
images  were  demonstrated. 
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The  signal -to-nol se  ratio  (SNR)  of  the  simulated  correlography 
data  was  measured  and  compared  with  theoretical  predictions,  and  It  was 
found  that  they  generally  agreed  well  with  one  another.  The  number  of 
frames  of  data  required  to  achieve  a  given  resolution  for  a  particular 
Image  of  a  satellite  model  was  determined.  However,  this  result 
depends  on  the  spatial -frequency  content  of  the  object.  The  effect  of 
photon  noise  for  low  light  levels  was  analyzed,  and  It  was  found  that 
the  effects  of  photon  noise  are  small  as  long  as  the  number  of  photons 
per  speckle  per  frame  Is  much  greater  than  two.  Alternative  Wiener 
filters  for  Improving  the  SNR  were  derived,  and  an  iterative  filtering 
method  was  suggested. 

■if  ,  • 

As  described  In  Section  4,  derivations  of  the  results  of 
noncoherent  averaging  of  Images  and  of  the  correlography  quantities  for 

'  *  i  i 

the  case  of  a  mixed  object  were  performed.  By  a  mixed  object  we  mean  a 
coherently  Illuminated  object  that  has  both  a  fixed,  deterministic 
component  (such  as  a  glint  or  glints)  and  a  random,  diffuse  component. 
It  was  found  that  If  the  deterministic  component  of  the  object  consists 
of  a  single  glint  (a  fairly  common  occurance),  then  the  traditional 
correlography  estimators  give  an  Incorrect  answer?  however,  a  new 
estimator,  denoted  as  <1^2  -  7  +  I2>,  or  <1  jl2  -  Var(l)>,  gives  the 
correct  answer.  This  new  estimator  also  gives  the  correct  answer  when 
no  glint  is  present. 

As- -  described  In  Section  5,  methods  were  1  developed  for 
reconstructing  a  coherent  Image  from  a  single  frame  of  aperture-plane 
Intensity  data  when  the  object  has  one  or  more  glints.  The  most 
successful  method  Is  effective  for  the  most  difficult  case  --  when  the 
object  has  multiple  glints  not  spatially  separated  frbm  the  diffuse 
component.  It  consists  of  three  successive  algorithms:  (1)  an 
autocorrelation  tri-intersection  algorithm  that  determines  the  glint 
positions  and  values;  (2)  the  AF-synthesis  algorithm  that  produces  a 
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partially-reconstructed  Image;  and  (3)  the  Iterative  Fourier  transform 
algorithm,  which  completes  the  reconstruction  of  a  fine-resolution 
Image.  This  method  was  demonstrated  to  reconstruct  fine-resolution 
Images  from  computer-simulated  data.  The  effect  of  noise  was  computer- 
simulated,  and  the  sensitivity  of  the  method  to  photon  noise  and  to 
glint  strength  was  determined  by  computer  reconstruction  experiments. 

Other  coherent  reconstruction  approaches  for  the  case  of  objects 
having  glints  were  also  Investigated.  Reconstructions  were 
successfully  performed  using  only  the  Iterative  Fourier  transform 
algorithm  for  the  cases  of  one  and  two  glints  on  the  object.  A. 
recursive  reconstruction  algorithm  based  on  the  autocorrelation  of  the 
object  (which  can  be  computed  from  the  aperture-plane  Intensity  data) 
was  developed  for  the  case  of  a  single  glint  on  the  object.  The 
effects  of  large  glints  on  the  quantization  error  when  detecting  the 
Intensity  data  were  analyzed.  We  found  that  the  quantization  error 
could  be  greatly  reduced  by  having  an  automatic  gain  control  that  would 
scale  down  the  detected  Intensity  when  a  very  large  glint  would  appear. 
A  variable  zero-offset  was  found  to  be  useful  to  a  lesser  extent. 

As  described  In  Sections  6  and  7,  methods  were  developed  for 
reconstructing  a  coherent  Image  from  a  single  frame  of  aperture-plane 
Intensity  data  when  partial  Information  about  the  phase  of  the  optical 
field  Is  available.  These  methods  are  also  effective  when  we  have 
partial  phase  Information  for  the  case  of  Incoherent-Image 
reconstruction  as  well. 

As  described  In  Section  6,  a  new  variation  of  the  Iterative  Fourier 
transform  algorithm,  called  the  expanding  weighted  modulus  algorithm, 
was  developed.  It  can  be  applied  when  the  Fourier  phase  Is  known  well 
over  a  small  aperture,  but  Is  unknown  over  the  large,  full  aperture. 
It  Involves  Iterating  with  progressively  larger  weighting  functions 
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Imposed  on  the  Fourier*  modulus  data,  reconstructing  progressively 
finer-resolutlon  Images,  effectively  bootstrapping  from  the  phase  dver 
the  small  aperture  to  the  phase  over  the  large  aperture.  The  Image- 
domain  constraint  used  for  the  Iterative  transform  algorithm  Is  a 
support  (finite  extent)  constraint. 

A  second  case,  reconstruction  from  one  bit  of  phase  known  over  the 
entire  aperture,  was  also  Investigated.  When  the  object  Is  "causal" 
(l.e.,  Is  entirely  to  one  side  of  the  optical  axis),  then  the  Image  can 
be  reconstructed  easily  by  a  windowing  operation  and  the  Iterative 
transform  algorithm. 

A  third  case  Is  where  the  phase  Is  known  poorly  over  the  entire 
aperture,  or,  equivalently,  a  noisy  phase  exists  over  the  entire 
aperture.  This  would  be  the  case  If  we  use  any  other  Imaging  method 
that  results  In  phase  errors.  For  this  case,  another  variation  of  the 
Iterative  transform  algorithm,  called  the  phase  variance  algorithm,  was 
developed,  In  It,  the  given  Fourier  modulus,  which  Is  assumed  to  be 
measured  with  a  higher  SNR,  Is  reinforced  exactly,  but  the  given 
Fourier  phase  Is  reinforced  Inexactly,*  it  Is  allowed  to  wander  from  the 
measured  phase  In  accordance  with  the  standard  deviation  of  the  error 
of  the  phase.  Appropriate  data  were  simulated  and  reconstruction 
experiments  were  performed.  For  the  case  of  Incoherent  Images  (for 
which  one  has  a  nonnegativity  constraint),  the  phase  variance  algorithm 
converges  faster  to  the  solution  than  the  traditional  algorithm.  For 
the  case  of  coherent  Images,  the  phase  variance  algorithm  produced 
Images  of  significantly  better  quality  than  that  given  by  the  noisy 
measured  phase,  but  left  room  for  further  improvement.  Thus,  the  phase 
variance  algorithm  should  be  used  to  clean  up  Images  produced  by  other 
Imaging  methods.  Investigation  of  the  effect  of  photon  noise  on  the 
Fourier  modulus  data  was  performed  by  computer  simulation  and 
reconstruction  experiments.  It  was  shown  that  If  the  Fourier  modulus 
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data  are  sufficiently  noisy,  then  the  phase  variance  algorithm  no 
longer  Improves  the  Image.  A  further  algorithm  improvement  was 
suggested:  allow  both  the  modulus  and  phase  to  vary  from  their  measured 
values,  each  In  accordance  with  the  variance  of  the  noise  on  that  data. 
Then  In  the  areas  of  spatial -frequency  space  where  the  modulus  Is  less 
noisy  than  the  phase,  the  phase  Is  Improved  by  reinforcing  the  modulus; 
and  In  the  areas  where  the  phase  Is  less  noisy  than  the  modulus,  the 
modulus  Is  Improved  by  reinforcing  the  phase. 

As  described  In  Section  7,  another  method  of  correcting  phase 
errors  over  the  entire  aperture  was  developed.  Called  2-D  shear 
averaging,  It  corrects  phase  errors  using  a  priori  Information  about 
the  statistics  of  the  coherent  object  --  that  It  Is  delta-correlated  — 
rather  than  the  usual  object-domain  support  constraint.  The  algorithm 
Is  computationally  simple,  consisting  of  three  steps.  First,  phase- 
error  differences  In  each  of  the  two  dimensions  are  estimated  from  the 
summation  over  a  sheared  product  of  the  given  Fourier  transform.  The 
derivation  of  this  result  Is  similar  to  that  of  the  shearing- 
interferometer  wavefront  sensor.  Next,  the  phase  error  Is  computed 
from  the  phase-error  differences  by  a  recruslve  method,  such  as  those 
employed  for  wavefront  sensing  or  Knox-Thompson  Image  reconstruction. 
Complex  exponentials  are  employed  to  avoid  phase  unwrapping 
difficulties.  Finally,  the  phase-error  estimate  Is  subtracted  from  the 
given  Fourier  phase  to  yield  a  corrected  Fourier  transform.  Inverse 
transformation  yields  the  corrected  Image.  Several  versions  of  the 
algorithm,  differing  In  the  way  that  the  sheared  product  Is  averaged, 
were  studied.  Analysis  was  performed  to  predict  the  residual  phase 
errors  left  by  this  statistically-based  algorithm.  Data  with  a  variety 
of  phase  errors  were  digitally  simulated  and  Image  reconstruction 
experiments  were  performed.  For  smooth  but  higher-order  2-D  phase 
errors  the  2-D  shear  averaging  algorithm  was  shown  to  Improve  the  Image 
quality  substantially,  although  there  was  room  for  further  Iniage- 
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quality  Improvement.  Thus,  the  2-D  shear  averaging  algorithm  should  be 
used  to  clean  up  Images  produced  by  other  Imaging  methods  when  the 
phase  errors  are  slowly  varying. 

As  described  In  Section  8,  several  different  classes  of  Imaging 
systems  appropriate  for  Imaging  satellites  were  briefly  compared, 
including  both  active  and  passive  approaches.  Several  methods  for 
obtaining  f Irresolution  Images,  besides  the  conventional  method*  of 
compensated  U*g<.jg  and  microwave  synthetic-aperture  radar,  appear  to 
show  promise.  One  good  example  of  a  novel  approach  Is  an  array  of 
sensors,  each  consisting  of  a  dual-plane  detection  of  laser  radiation 
plus  a  wavefront  sensor.  This  approach  employs  a  robust  form  of  phase 
retrieval  to  determine  the  optical  field  In  the  aperture  plane  from  two 
arrays  of  Intensity  data  without  requiring  heterodyne  detection.  The 
wavefront  sensor,  operating  on  Incoherent  light,  determines  the 

'  wavefront  error  due  to  atmospheric  turbulence.  The  phase  of  the  latter 
Is  subtracted  from  the  phase  of  the  former  to  yield  the  phase  due  to 
the  object  alone. 

In  summary,  several  different  approaches  to  reconstructing  fine- 
resolution  Images  of  satellites,  using  relatively  simple  receivers, 
were  developed.  These  methods  could  make  possible  Imaging  systems  of 
greatly  reduced  cost  and  complexity  compared  with  compensated  Imaging 
(using  adaptive  optics).  They  tend  to  scale  well  for  the  case  of  deep- 
space  objects,  for  which  the  receiver  array  must  be  much  larger  than 
any  existing  optical  telescope. 
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3.0  IMAGING  CORRELOGRAPHY 


In  this  section,  the  new  Imaging  method  we  call  Imaging 
correlography  Is  described.  It  makes  use  of  multiple  realizations  of 
the  Intensity  of  the  coherent  optical  field,  backscattered  from  the 
object,  measured  In  the  aperture  plane,  to  arrive  at  an  Incoherent 
Image  of  the  object.  The  basic  concept  Is  described  In  Section  3.1. 
The  method  Is  demonstrated  for  sparse  arrays  of  detectors  In  Section 
3.2.  The  slgnal-to-nolse  ratio  achieved  In  computer  simulation  Is 
compared  with  theory  In  Section  3.3.  Wiener  filtering  Issues  are 
further  discussed  In  Section  3.4.  Correlography  for  a  "mixed"  object, 
l.e.  one  having  a  deterministic  (glint)  component  as  well  as  a  diffuse 
component,  Is  described  In  Section  3.5,  where  It  Is  shown  that  a  new 
estimator  Is  required. 

3.1  IMAGING  CORRELOGRAPHY  THEORY  AND  RESULTS 

It  Is  well  known  that  the  spatial  structure  of  a  fully  developed 
laser-speckle  pattern  —  produced  by  the  coherent  Interference  of  light 
backscattered  from  a  sufficiently  diffuse,  reflecting  surface  —  Is 
dependent  on  the  macroscopic  features  of  the  Illuminated  surface  [3.1]. 
In  this  Chapter  we  demonstrate  that  measurements  of  the  backscattered 
speckle  Intensity  are  sufficient  to  (uniquely)  reconstruct  a  high- 
resolution,  unspeckled,  Incoherent  Image  (or  brightness  distribution) 
of  the  coherently  Illuminated  object. 

Our  approach  to  Image  synthesis  Is  based  on  the  fact  that  from  the 
average  energy  spectrum  of  a  laser-speckle  Intensity  pattern  we  can 
obtain  the  autocorrelation  function  of  the  Illuminated  object's 
brightness  distribution  [3.2].  Here,  the  object's  brightness 
distribution  corresponds  to  the  object's  reflectance  function  or, 
alternatively,  to  Its  Irradlance  distribution  had  the  object  been 
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Illuminated  by  an  Incoherent  light  source.  Since  the  Fourier  transform 
of  the  autocorrelation  of  the  object  brightness  function  Is  equivalent 
to  the  squared  modulus  of  the  Fourier  transform  of  the  brightness 
function  [3.3],  an  Image  of  the  object  can  be  obtained  If  the  phase 
associated  with  this  Fourier  transform  can  be  determined,  Fortunately, 
a  practical  solution  to  this  phase-retrieval  problem  has  been 
demonstrated  by  Flenup  [3. 4-3.6],  In  which  an  Iterative  transform 
algorithm  can  be  used  to  recover  the  phase  associated  with  the  modulus 
of  the  Fourier  transform  of  a  real,  nonnegative  object  function, 
provided  that  certain  boundedness  and  nonnegativity  constraints  are 
continually  reinforced  throughout  the  Iteration  process.  The  Iterative 
transform  algorithm,  together  with  certain  digital  preprocessing 
operations  (which  are  described  below)  permit  us  to  recover  complete, 
unspeckled  Images  from  non Imaged  speckle  data. 

Let  us  suppose  that  a  diffuse  object  Is  flood  Illuminated  with  a 
laser  whose  coherence  length  Is  at  least  twice  as  long  as  the  object  Is 
deep.  An  array  of  photodetectors  measures  the  backscattered  light 
Intensity  In  a  far-fleld  plane  some  distance  from  the  object.  We 
assume  that  the  object  Is  optically  rough,  so  that  Its  microscale 
surface  height  variations  are  random  and  comparable  In  size  with  the 
wavelength  of  light.  This  being  the  case,  the  reflected  laser  light  Is 
randomly  (and  coherently)  dephased,  and  the  photodetectors  In  the 
observation  plane  record  a  fully  developed  laser-speckle  pattern. 

Each  realization  of  the  observed  speckle  Intensity  In(u)  may  be 
expressed  as  the  squared  modulus  of  the  Fourier  transform  of  the 
complex  object  field: 

In(u)  -  IFn(u)12  -  iaF</n(*)>l2  ,  (3-1) 
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where  7  denotes  a  Fourier  transform,  /n(x)  ■  l/Q(x)l  exp[i#n3 ~  is  the 
field  reflected  by  the  object,  l/0(x)l  is  the  object's  field  amplitude 
reflectivity,  and  #n(x)  is  the  (random)  phase  of  the  nth  realization  of 
the  reflected  object  field  associated  with  the  object's  surface  height 
profile.  In  the  above  expression,  x  represents  a  two-dimensional 
spatial  (or  angular)  coordinate  vector  in  object  or  image  space;  u 
represents  a  two-dimensional  coordinate  in  the  measurement  plane.  The 
Inverse  Fourier  transform  of  the  observed  speckle  pattern  is 
proportional  to  the  autocorrelation  of  the  object  field,  which  may  be 
written  as 


rn(x)  -  7{IFn(u)l  H(u)> 

-  [/„(*)  •  /„(*)]  *  h(x)  ,  (3-2) 

where  7"1  denotes  an  Inverse  Fourier  transform,  *  denotes  a  convolution 
operation,  and  9  denotes  an  autocorrelation.  The  aperture  function 
H(u)  denotes  the  region  of  the  measurement  plane  over  which  the  speckle 
pattern  Is  observed:  H(u)  ■  1  for  points  within  the  measurement 
aperture,  and  H(u)  ■  0  elsewhere.  The  function  h(x)  ■  T1{H(u)}  is  the 
(diffraction-limited)  coherent  impulse  response;  hence  rn(x)  is  a 
diffraction-limited  (albeit  speckled)  autocorrelation  of  the  laser- 
illuminated  object. 

Using  the  Iterative  transform  algorithm,  one  could  attempt  to 
reconstruct  a  complex-valued,  speckled  image  of  /n(x)  from  1 Fn (u) 1 2 
H(u)  or  equivalently  from  rn(x).  However,  at  present  the  practical 
reconstruction  algorithm  is  effective  only  for  certain  classes  of 
complex-valued  objects  if  the  object's  support  is  known  a  priori  [3.7] 
and  for  even  more  restrictive  classes  of  complex-valued  objects  If  the 
object  support  is  unknown.  (The  support  is  the  closed  set  of  points 
outside  which  the  object  is  zero.)  Such  cases  are  described  in  Section 
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4.  In  this  section,  we  concentrate  on  a  method  that  allows  us  to 
reconstruct  a  real,  nonnegative  image  —  a  case  for  which  the  iterative 
transform  algorithm  is  effective  for  a  broad  class  of  objects. 

Image  recovery  begins  by  estimating  the  average  energy  spectrum  of 
the  observed  speckle  pattern  by  averaging  together  the  squared  moduli 
of  many  Independent  speckled  autocorrelations  rp(*).  This  may  be 
referred  to  as  noncoherent  averaging  of  the  coherent  autocorrelations. 
Independent  realizations  of  the  speckle  pattern  can  be  obtained,  for 
example,  by  laterally  displacing  the  observation  plane  with  respect  to 
the  object  or  by  measuring  the  speckle  pattern  for  slightly  different 
rotations  of  the  object.  We  can  show  that  as  the  number  N  of 
Independent  speckle  realizations  Increases,  the  computed  average  energy 
spectrum  converges  to  [3.8] 

11m  N"1  23  !rnW2  "  blh(x)!2  +  cr0(x)  *  I h(x) 1 2  ,  (3-3) 

N-m*  n^l  n  0 

where 

2 

b  ■  c[  J  l/0(x')l2  d2x'j  (3-4) 

Is  the  square  of  the  total  measured  irradiance, 

f,(»l  ■  l/0W|!  •  »-5) 

Is  the  autocorrelation  of  the  object  brightness  function,  and  c  Is  a 
constant.  Thus  the  average  energy  spectrum  converges  to  the  sum  of  an 
autocorrelation  of  the  desired  Incoherent  image  plus  a  dc  term 
b I h (x) 1 2 ,  where  the  dc  term  Is  simply  the  (Incoherent)  polntspread 
function  of  the  collecting  aperture,  possessing  a  strength  b.  On 
subtracting  the  dc  term  from  the  averaged  energy  spectrum,  we  obtain  a 
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diffraction-limited  autocorrelation  of  the  Incoherent  object.  The 
square  root  of  the  Fourier  transform  of  this  incoherent  object 
autocorrelation,  then,  provides  us  with  an  estimate  of  the  modulus  of 
the  Fourier  transform  of  the  object's  brightness  function.  Note  that 
one  can  obtain  the  same  results  by  subtracting  a  bias  from  an  average 
of  the  autocorrelations  of  In(u)  and  then  taking  the  square  root.  One 
can  see  that  the  latter  approach  Is  analogous  to  a  highly  redundant, 
multichannel  Intensity  Interferometer  [3.9].  This  latter  approach  Is 
described  In  more  detail  In  Section  3.2. 

We  conducted  a  series  of  computer  experiments  to  demonstrate  that 
phase  retrieval  can  be  used  to  recover  Imagery  from  speckle  data 
processed  In  this  way.  Original  object  data  for  these  experiments  were 
contained  In  o  digitized  photograph  of  a  satellite  model  Illuminated 
with  Incoherent  light.  These  data  comprised  approximately  40  x  60 
pixels  embedded  In  a  128  x  128  discrete  array.  Each  realization  of  a 
coherent  (speckled)  Image  of  the  object  was  obtained  from  the  digitized 
photograph  by  (1)  replacing  each  pixel  with  a  circular-complex  Gaussian 
random  variable  whose  real  and  Imaginary  parts  possessed  variances 
equal  to  half  of  the  pixel  Intensity  value  and  (2)  low-pass  filtering 
the  result.  The  filter  used  to  smooth  the  complex  object  data 
corresponds  to  the  aperture  function  H(u),  which  was  represented  by  a 
64  x  64.  square  of  detector  pixels  embedded  In  128  x  128  array. 
Multiple  realizations  of  the  coherent  object  data  were  obtained  by 
using  different  random-number  seeds  In  the  computation  of  the  complex 
Gaussian  random  variables.  Each  coherent  Image  autocorrelation  rn(x) 
was  computed  by  Inverse  Fourier  transforming  the  squared  modulus  of  the 
apertured  Fourier  transform  of  the  simulated  coherent  Image.  Averages 
of  both  the  speckled  autocorrelations  and  their  squared  moduli  (l.e., 
the  energy  spectrum  of  the  speckle-intensity  patterns)  were  then  taken. 
A  function  proportional  to  the  square  of  the  former,  an  estimate  of  the 
dc  term,  was  subtracted  from  the  latter  (the  noncoherent  average)  to 
arrive  at  an  estimate  for  the  autocorrelation  of  the  Incoherent  Image. 
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The  process  of  noncoherently  averaging  object-field 
autocorrelations  and  subtracting  the  dc  term  is  illustrated  in  rig. 
3-1.  The  first  column  contains  averages  of  the  squared  inverse  Fourier 
transforms  of  N  simulated  speckle  measurements  providing  estimates  of 
the  speckle  energy  spectrum,  where  N  is  the  number  of  independent 
speckled  autocorrelations  noncoherently  averaged.  The  second  column 
shows  the  corresponding  dc  term,  which,  for  the  case  of  a  square 
aperture,  is  a  squared  sinc(x)  [i.e.,  (rx)"*  sin(rx)]  function.  The 
third  column  shows  the  results  when  the  dc  term  is  subtracted  from  the 
noncoherently  averaged  autocorrelations  of  the  first  column.  Note  that 
the  speckle  artifacts  in  the  averaged  autocorrelations  (in  the  first 
and  third  columns  of  Fig.  3-1)  disappear  as  N  increases. 

The  incoherent  autocorrelation  estimate  (with  the  dc  term  removed) 
was  then  Fourier  transformed,  and  the  square  root  was  taken,  to  arrive 
at  an  estimate  of  the  modulus  of  the  Fourier  transform  of  the  object 
brightness  function.  Negative  numbers,  resulting  from  noise  associated 
with  the  finite-average  approximation  to  an  ensemble  average,  were  set 
to  zero  before  the  square  root  was  taken.  Images  were  reconstructed 
from  the  Fourier  modulus  estimates  by  using  the  Iterative  Fourier- 
transform  algorithm  [3.5,  3.6],  using  several  cycles  of  the  hybrid 
input-output  algorithm  (using  p  -  0.7)  and  the  error  reduction 
algorithm  until  the  algorithm  appeared  to  stagnate.  The  object-domain 
constraints  used  were  nonnegativity  (since  an  incoherent  image  is  being 
reconstructed)  and  a  loose  support  constraint  (a  rectangle  half  the 
size  of  the  smallest  rectangle  enclosing  the  autocorrelation). 

Data  along  the  first  row  of  Fig.  3-2  Illustrate  a  direct 
application  of  the  phase-retrieval  algorithm  to  the  Fourier  modulus 
estimate.  Figure  3-2(A)  represents  the  dc-subtracted  autocorrelation 
for  N  ■  10^  Independent  speckle  patterns.  Figure  3-2 (B)  shows  the 
corresponding  Fourier  modulus  data  produced  by  Fourier  transforming  the 
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FIGURE  3-1.  ESTIMATING  THE  ENERGY  SPECTRUM  OF  SPECKLE  INTENSITY  BY 
NONCOHERENTLY  AVERAGING  MANY  COHERENT  SPECKLED  IMAGE  AUTOCORRELATIONS. 
(A)  Noncoherent  average  of  N  =  4  autocorrelations;  (B)  estimate  of  dc 
term;  (C),  (A)  minus  (B);  (D)-(F)  N  =  32;  (G)-(I)  N  *  128;  (J)-(L)  N  - 
1024. 
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FIGURE  3-2.  IMAGE  RECOVERY  FROM  NONCOHERENTLY  AVERAGE  AUTOCORRELATION 
DATA  (N  -  10,000,  FILLED  APERTURE).  (A)  DC-adjusted,  noncoherently 
averaged  autocorrelations,  (B)  estimate  of  the  Fourier  modulus  of  the 
incoherent  object,  (C)  image  reconstructed  from  (B)  using  the  iterative 
transform  (phase-retrieval)  algorithm,  (D)  Wiener  filter,  (E)  filtered 
Fourier  modulus  estimate,  (F)  image  reconstructed  from  (E),  (G) 
original  incoherent  object,  (H)  Wiener  filtered,  incoherent  object,  (l) 
result  of  Wiener  filtering  (C). 


averaged  autocorrelation  [Fig.  3-2 (A)]  and  then  taking  the  square  root. 
Figure  3-2(C)  Is  the  reconstructed  Image  produced  by  applying  the 
phase-retrieval  algorithm  as  outlined  above.  Note  that  this  Image  is 
very  noisy  compared  with  the  original  Incoherent  object,  shown  In  Fig. 
3-2(G).  Noise  In  the  reconstructed  Image  Is  due  to  the  fact  that  a 
finite  number  of  speckle  realizations  were  used  to  estimate  the  Fourier 
modulus.  To  reduce  these  noise  effects,  we  multiplied  the  Fourier 
modulus  estimate  [Fig.  3-2(B)]  by  a  Wiener  filter  of  the  form 


,  v  OTF(Au)  E,(Au) 

W(Au)  «  - - w--5 -  , 

I0TF(A) «  Es(Au)  +  En 


(3-6) 


where  OTF(Au)  -  H(u)9H(u)  Is  the  optical  transfer  function  of  the 
receiver  aperture,  E$(Au)  Is  an  average  energy  spectrum  for  objects  of 
this  type  (estimated  by  taking  an  angular  average  over  the  squared 
Fourier  modulus  of  the  object),  and  En  Is  the  energy  spectrum  of  the 
noise.  We  approximated  En  by  a  constant  whose  value  was  obtained  by 
averaging  the  squared  Fourier  modulus  estimate  over  those  higher 
spatial  frequencies  where  the  slgnal-to-noise  ratio  was  less  than  one. 
Figure  3-2(D)  shows  the  Wiener  filter  used  for  this  example. 

Figure  3-2(E)  shows  the  filtered  Fourier  modulus  estimate  equal  to 
the  product  of  Figs.  2(B)  and  2(D).  Figure  3-2(F)  shows  the  Image 
reconstructed  from  the  Wiener-filtered  Fourier  modulus  estimate  using 
the  phase-retrieval  algorithm.  Note  that  the  Wiener  filter  has 
significantly  Improved  the  quality  of  the  reconstructed  Image  In  Fig. 
3-2(F)  over  that  In  Fig.  3-2(C)  reconstructed  without  Wiener  filtering. 
For  comparison,  the  original  object  [shown  In  Fig.  3  2(G)]  was  passed 
through  the  Wiener  filter  of  Fig,  3«2(D),  with  the  result  shown  In  Fig. 
3-2(H).  The  Image  reconstructed  from  speck le-correl awl  on  measurements, 
shown  In  Fig.  3-2(F),  compares  favorably  with  the  filtered  object  [Fig. 
3-2{H) ] ,  Indicating  good  performance  on  the  part  of  the  Iterative 
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transform  reconstruction  algorithm.  Finally,  Fig.  3-2(1)  shows  the 
result  of  applying  the  Wiener  filter  to  the  reconstructed  Image  shown 
in  Fig.  3-2(C).  Apparently,  Wiener  filtering  followed  by  Image 
reconstruction  is  superior  to  image  reconstruction  followed  by  Wiener 
filtering. 

» 

These  results  demonstrate  the  possibility  of  recovering  Images  from 
nonlmaged  laser  speckle  patterns:  by  averaging  over  many  realizations 
of  the  coherent  (speckle)  Intensity  data,  an  estimate  of  the 
autocorrelation  and  Fourier  modulus  of  the  Incoherent  object  can  be 
obtained.  And,  from  the  Fourier  modulus  estimate.  It  Is  possible  to 
reconstruct  an  unspeckled  Image  by  applying  a  phase-retrieval  algorithm 
with  a  nonnegativity  constraint. 

Figures  3-3  and  3-4  show  results  similar  to  those  in  Fig.  3-2,  ,but 
for  a  smaller  number  of  realizations,  N,  of  the speckle  patterns. 
Figure  3-3  shows  the  case  for  N  ■  1024  realizations  and  Figure  3-4  the 
case  of  N  ■  128  realizations.  As  expected,  the  Image  quality  decreases 
with  fewer  number  of  speckle  frames  due  to  the- statistical  noise 
associated  with  a  finite  number  of  frames,  Section  3-3  will  discuss  In 
greater  detail  the  slgnal-to-nolse  Issues. 

Since  these  computer  simulations  were  performed,  laboratory 
experimental  verification  of  Imaging  correlography  has  also  been 
accomplished  [3.10]. 

This  section  Is  an  expansion  of  Reference  3.11. 
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FIGURE  3-3.  IMAGE  RECOVERY  FROM  NONCOHERENTLY  AVERAGE  AUTOCORRELATION 
DATA  (N  -  1024,  FILLED  APERTURE).  (A)  DC-adjusted,  noncoherently 
averaged  autocorrelations,  (B)  estimate  of  the  Fourier  modulus  of  the 
incoherent  object,  (C)  image  reconstructed  from  (B)  using  the  iterative 
transform  (phase-retrieval)  algorithm,  (D)  Wiener  filter,  (E)  filtered 
Fourier  modulus  estimate,  (F)  image  reconstructed  from  (E),  (G) 
original  incoherent  object,  (H)  Wiener  filtered,  incoherent  object,  (I 
result  of  Wiener  filtering  (C). 


FIGURE  3-4.  IMAGE  RECOVERY  FROM  NONCOHERENTLY  AVERAGE  AUTOCORRELATION 
DATA  (N  ■  128,  FILLED  APERTURE).  (A)  DC-adjusted,  noncoherently 
averaged  autocorrelations,  (B)  estimate  of  the  Fourier  modulus  of  the 
incoherent  object,  (C)  Wiener  filter,  D)  filtered  Fourier  modulus 
estimate,  (E)  image  reconstructed  from  (D),  (F)  original  incoherent 
object,  (G)  Wiener  filtered,  incoherent  object. 
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3.2  IMAGING  CORRELOGRAPHY  WITH  SPARSE  ARRAYS  OF  DETECTORS 

In  this  section  the  use  of  Imaging  correlography,  Introduced  In  the 
previous  section,  with  sparse  arrays  of  detectors  is  discussed  and 
demonstrated  through  digital  simulations.  In  this  case  It  Is  important 
to  emphasize  the  relationship  between  the  aperture  function  shape  and 
the  modulation  transfer  function  (MTF)  for  the  Image.  For  this  reason 
we  start  with  an  alternative  (but  mathematically  equivalent) 
explanation  of  the  correlography  process. 

Rather  than  relating,  as  we  have  done  above  In  Eq.  (3-3),  the 
average  energy  spectrum  of  the  speckle  pattern  to  the  autocorrelation 
function  of  the  object's  brightness  function,  we  can  equate  the 
autocovariance  of  the  far-fleld  laser  speckle  pattern  with  the  energy 
spectrum  of  the  object's  brightness  function.  This  second 
interpretation  suggests  the  following  procedure  for  Image  recovery: 
(1)  estimate  the  autocovariance  of  the  observed  speckle  Intensity,  (2) 
take  the  square  root  of  the  estimated  autocovariance,  (3)  recover  the 
phase  associated  with  this  square-root,  and  finally  (4)  Inverse  Fourier 
transform  the  assembled  Fourier  data.  Image  recovery  using  this 
prescription  uncovers  a  close  relationship  between  Imaging 
correlography  and  Image  recovery  from  Intensity  Interferometry  [3.13], 
where  the  object's  1  Fourier  phase  Information,  too,  Is  lost  to  the 
measurement  process.  (The  fact  that  Fourier  domain  Information  of 
Incoherent  objects  can  be  obtained  from  far-fleld  correlations  Is,  of 
course,  a  consequence  of  the  Van  Cittert-Zernlke  theorem  [3.12].) 

We  can  demonstrate  the  relationship  between  the  autocovariance  of 
the  laser  speckle  pattern  and  the  object  energy  spectrum  by  considering 
the  measurement  process  Involved  In  Imaging  correlography.  Let  us 
suppose  that  a  diffuse  object  Is  flood  Illuminated  with  a  laser  so  that 
the  object  lies  entirely  within  the  coherence  volume  of  the  laser  beam. 
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A  two-dimensional  array  of  photodatectors  measures  the  backscattered 
light  intensity  In  a  (far-fleld)  plane  some  distance  z  from  the  object 
(see  Fig.  3-5  for  a  possible  measurement  scenario).  We  assume  that  the 
object  Is  optically  rough  so  that  Its  microscale  surface  height 
variations  are  random  and  of  size  comparable  to  or  greater  than  the 
wavelength  of  light.  Additionally,  we  assume  that  the  transverse  scale 
size  of  the  surface  height  fluctuations  Is  small  compared  to  the 
resolution  patch  size  associated  with  the  collecting  array  (l.e.,  the 
spatial  correlation  of  surface  roughness  is  small  compared  to  Xz/D, 
where  X  Is  the  wavelength  of  light,  z  Is  the  range,  and  D  Is  the 
largest  array  dimension).  This  being  the  case,  the  reflected  laser 
light  Is  randomly  (and  coherently)  dephased  and  the  photodetectors  In 
the  observation  plane  record  the  Intensity  pattern  of  a  fully  developed 
laser  speckle  pattern  [3.1] *  . 

An  estimate  of  the  autocovariance  of  the  measured  speckle  pattern, 
In  [see  Eq .  3-1] ,  may  be  computed  as  follows  from  N  realizations  of  the 
laser  speckle  Intensity: 

Cj(Au;  N)  i  |  ^  JJ  H(u  +  Au)  H(u)  [lR(u  +  Au)  ln(u)  -  T2]  d2u 

-  jj  H(u  t  Au)  H(u)  |j|  £  [ln(u  +  Au)  In(u)  -  T2]}  d2u  (3-7) 

where  T  Is  the  average  Intensity  of  the  observed  speckle  pattern,  Au  Is 
a  vector  separation  In  the  measurement  plane,  and  H(u)  Is  the  pupil 
(aperture)  function  associated  with  the  collecting  array,  defined  as 


H(u) 


1, 

.0, 


for  u  e  aperture  array 
otherwise 


(3-8) 
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In  the  limit  as  N  (the  number  of  Independent  observed  speckle  patterns) 
approaches  Infinity,  one  can  use  the  moment  factoring  theorem  for 
circular-complex  Gaussian  (ccg)  fields  [3.13]  to  show  that 

■iJ?*  I  j{j!j  [!n(u  +  Au)  !n<u)  '  ?2]  “  irUu)l2  (3-9) 

where  r(Au)  ■  ^[1 fQ(x) I  ]  Is  the  Fourier  transform  of  the  object's 
brightness  distribution  [1.e.(  r(Au)  Is  the  mutual  intensity  of  the 
(complex)  speckle  field  in  the  measurement  aperture,  evaluated  at  field 
points  separated  by  a  vector  Au].  Our  ability  to  Invoke  the  ccg  moment 
theorem  above  follows  from  the  fact  that  the  observed  speckle  field  Is 
ccg,  since  the  speckle  pattern  Is  fully  developed.  In  the  limit  N  ♦  •», 
we  therefore  find  from  Eqs.  (3-7)  and  (3-9)  that  the  estimated 
autocovariance  of  the  speckle  Intensity  observed  over  the  measurement 
aperture  H(u)  Is  given  by 

Ct(Au)  i  11m  Cr(Au;  N) 

1  N-m»  1 

-  OTF(Au)  ir(Au)l2  (3-10) 

where  OTF(u)  Is  the  autocorrelation  of  H(u).  This  result  demonstrates 
that  Cj(Au;N)  provides  an  estimate  for  ir(Au)l2,  the  energy  spectrum  of 
the  object's  brightness  function  —  the  square  root  of  which  Is  an 
estimate  of  the  Fourier  modulus  of  the  object's  brightness  function. 
This  square  root  Is  used  In  the  Iterative  transform  algorithm  to 
retrieve  the  associated  Fourier  phase  data  and  thereby  reconstruct  an 
Image. 
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We  see  from  Eq.  (3-10)  that  the  estimated  autocovariance  of  the 
observed  speckle  pattern  provides  a  weighted,  or  filtered,  estimate  of 
the  object's  energy  spectrum.  This  weighting  Is  completely  determined 
by  the  spatial  arrangement  of  the  detectors  making  up  the  collecting 
aperture.  Because  this  weighting  function  OTF(Au)  Is  equal  to  the 
autocorrelation  of  the  measurement  pupil,  we  refer  to  OTF(Au)  as  the 
optical  transfer  function  (OTF)  for  the  Imaging  correlography  system  — 
with  obvious  analogy  to  the  OTF  arising  In  the  analysis  of  Incoherent 
Imaging  systems.  The  modulation  transfer  function  (MTF),  Is  just  the 
modulus  of  the  OTF}  and  since  the  OTF  Is  nonnegative,  MTF(Au)  ■ 
OTF(Au).  The  fact  that  this  OTF  Is  In  the  form  of  an  autocorrelation 
allows  us  to  consider  the  use  of  sparse  arrays  of  Intensity  detectors 
In  Imaging  correlography. 

The  fact  that  the  OTF  for  Imaging  correlography  Is  given  by  the 
autocorrelation  of  the  pupil  function  H(u)  suggests  a  procedure  with 
which  to  remove  sldelobe  artifacts  Introduced  by  a  multiple-aperture 
(sparse  array)  measurement  scheme.  If  the  detector  elements  are 
positioned  so  that  the  autocorrelation  of  the  detector  array  does  not 
drop  to  zero  within  the  bandpass  of  the  OTF,  the  object  energy  spectrum 
estimated  by  the  Imaging  correlography  process  contains  essentially  the 
same  spatial  frequencies  as  a  filled  aperture  having  the  same  diameter 
as  the  sparse  array.  And,  provided  that  the  noise  In  the  estimated 
autocovariance  Is  not  too  great,  the  energy  spectrum  estimate  can  be 
boosted  to  match  the  OTF  of  a  completely  filled  aperture;  an  Image  with 
nearly  the  resolution  of  the  full  aperture  Is  then,  In  theory, 
synthesized. 

In  practical  applications  of  Imaging  correlography,  noise  In  the 
Fourier  modulus  estimate  will  arise  from  many  sources  Including 
detector  noise,  background  flux  noise,  photon  shot  noise,  and  noise 
that  Is  Introduced  when  a  finite  number  of  sreckle  measurements  Is  used 
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to  estimate  the  speckle  autocovariance.  Where  all  the  noise  sources 
additive  and  uncorrelated  with  the  signal  component,  one  would 
logically  Implement  the  MTF  boosting  procedure  by  applying  a  Wlener- 
Helstrom  filter  [3.14]  to  the  Fourier  modulus  estimate  so  that  the 
mean-square  error  between  the  estimated  Image  and  the  true  (full- 
resolutlon)  Image  Is  minimized.  Even  If  the  signal  and  noise  sources 
do  not  exactly  satisfy  these  conditions,  a  Wlener-Helstrom,  filter  Is 
still  very  advantageous  to  use  t3.14]. 


For  conventional  Incoherent  Imaging  systems  the  Wlener-Helstrom 
filter  Is  of  the  form 


W(4u) _ PTFltM)  irttjj)|2 _ _ 

IOTF(4u) lz  ir(4u)l'  +  E„(4u) 


(3-11) 


2 

where  I r (Au) I  Is  the  energy  spectrum  of  the  object's  brightness 
function,  equal  to  E$(Au)  used  In  Eq.  (3-6),  OTF(Au)  Is  the  OTF  of  the 
collecting  aperture,,  and  En(Au)  Is  the  energy  spectrum  of  the  Image- 
domain  noise.  This  filter  Is  based  on  a  model  of  the  Imaging  process 
which  Is  given  In  the  Fourier  domain  as  OTF(Au)  T(Au)  +  r'o1se. 
However,  a  better  model  for  Imaging  correlography  Is 


Cj (Au;  N)  -  OTF(Au)  IT(Au) I2  +  Nc(Au)  ,  (3-12) 


where  Nc(Au)  Is  additive  noise,  for  which  the  appropriate  filtering 
operation  to  estimate  the  object's  energy  spectrum  Is 

if(Au)!2  -  Wc(Au)  Cj(Au;  N)  (3-13) 
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where  the  filter  Is  given  by 


W  (iu) _ OTF(tu)  ir(tu)l4 

IOTF(Au)lz  ir(Au)l4  +  Er(Au) 

where  Ec(Au)  Is  the  variance  of  Nc(Au). 


(3-14) 


Whether  taking  the  square  root  of  the  speckle  autocovariance  then 
filtering  with  Eq.  (3-11)  or  filtering  the  speckle  autocovariance  with 
Eq.  (3-14),  then  taking  the  square  root,  In  either  case  the  MTF  Is 
boosted  where  the  slgnal-to-nolse  ratio  Is  high  and  It  Is  depressed 
where  the  slgnal-to-nolse  ratio  Is  low,  thereby  resulting  In  a  better 
Fourier  modulus  estimate.  Indeed,  results  of  the  computer  simulations 
presented  In  the  next  section  demonstrate  that  such  filtering 
techniques  Improve  the  overall  quality  of  Imagery  recovered  In  Imaging 
correlography. 


We  conducted  a  series  of  computer  experiments  to  demonstrate  that 
phase  retrieval  can  be  used  to  recover  Imagery  from  far-fleld  speckle 
Intensity  data  collected  over  a  sparse  array.  The  procedure  followed 
here  Is  essentially  the  same  as  that  reported  In  Section  3.1,  with  the 
exception  that  here  the  speckle  realizations  used  to  compute  an 
estimate  of  the  Incoherent  object's  energy  spectrum  are  masked  with  a 
pupil  function  H(u)  emulating  a  sparse  collecting  array. 

The  original  object  data  for  these  experiments  was  the  same  as 
described  In  Section  3.1.  For  these  sparse  aperture  simulations,  we 
used  a  Golay-type  array  [3.15]  comprising  six  subapertures.  Figure  3-6 
shows  the  Golay  aperture  configurations  used  for  this  study  together 
with  the  corresponding  OTF's  arid  point-spread  functions.  The  narrower- 
and  wlder-segment  Golay  arrays  were  both  configured  to  have  a  16  pixel 
separation  between  adjacent  suoapertures;  the  diameters  of  the 
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FIGURE  3-5.  GOLAY  CONFIGURATIONS  CONTAINING  SIX  SUBAPERTURES.  Upper 
left:  Aperture  functions  H(u)  for  the  Golay-6.  Right:  OTF's 
corresponding  to  the  Golay-6  aperture  functions  shown  in  the  upper 
left.  Lower  left:  Point-spread  function  associated  with  the  wider- 
segment  Golay-6  aperture. 
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individual  subapertures  In  the  narrower-  and  wider-segment  arrays  were 
11  and  13  pixels,  respectively.  In  both  the  cases  the  OTF,  which  Is 
the  autocorrelation  o*  the  pupil  function,  consists  of  a  large  central 
peak  surrounded  by  30  satellite  peaks.  Although  the  widths  of  the 
subapertures  for  both  cases  were  chosen  to  be  large  enough  that  the  OTF 
does  not  drop  to  zero  within  the  bandpass,  the  narrower-segment  array 
OTF  does  drop  to  low  values  in  the  regions  between  the  satellite  OTF 
peaks.  In  the  presence  of  noise,  these  dips  in  the  OTF  could  result  in 
infornation  loss  at  these  spatial  frequencies.  For  the  wider-segment 
case,  the  OTF  stays  above  half  of  the  value  of  the  satellite  peaks  In 
the  areas  between  the  satellite  peaks,  as  can  be  seen  In  Figure  3-7. 
For  this  reason  the  wider-segment  Golay  array  was  chosen  for  the 
simulation.  To  perform  the  filtering  operation  on  the  complex  object 
data,  the  sampled  Golay  arrays  were  embedded  in  a  128  x  128  array. 
Multiple  realizations  of  the  coherent  object  data  are  then  obtained  by 
using  different  random  number  seeds  in  the  computation  of  the  complex 
Gaussian  random  variables. 

An  estimate  of  the  object  energy  spectrum  was  formed  by  processing 
multiple  arrays  of  pupil-plane  speckle  intensity  data  computed  from 
realizations  of  the  filtered  coherent  object.  Several  different 
estimators  of  the  object  energy  spectrum  can  be  used,  such  as  the  one 
given  by  Eq.  (3-7).  Figure  3-8  shows  an  example  of  the  data  that  was 
computed  for  these  experiments,  first  the  average  energy  spectrum  of 
the  speckle  intensity  was  computed  by  Inverse  Fourier  transforming  the 
square  modulus  (i.e.,  the  speckle  intensity)  of  the  Golay-apertured 
Fourier  transform  for  each  simulated  coherent  image,  and  then  these 
speckled  energy  spectra  were  averaged  together  as  shown  in  Figure 
3-8 (A) .  After  a  large  number  N  of  independent  coherent  speckle  data 
sets  were  processed  in  this  fashion,  a  DC  term  (in  fact  a  function, 
corresponding  to  a  scaled  version  of  the  squared  modulus  of  an  average 
of  the  Fourier  transforms  of  the  windowed,  speckle  intensity  arrays 
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FIGURE  3-8.  ESTIMATING  THE  ENERGY  SPECTRUM  OF  SPECKLE  INTENSITY  FOR 
THE  GO LAY -6  APERTURE.  (A)  Noncoherent  average  of  1024  autocorrela¬ 
tions;  (B)  an  estimate  of  the  dc-term;  (C)  (A)  minus  (B);  (D)  the 
corresponding  estimate  of  the  Fourier  modulus  (the  square  root  of  the 
estimated  autocovariance). 


observed  over  the  measurement  aperture)  was  computed,  as  shown  In 
Figure  3-8(B).  This  DC  term  was  subtracted,,  with  the  result  shown  In 
Figure  3-8(C).  This  result  was  Fourier  transformed,  providing  an 
estimate  of  the  autocovariance  of  the  observed  speckle  pattern,  which 
by  Eq.  (3-10),  is  an  OTF-welghted  estimate  of  the  Incoherent  object's 
energy  spectrum.  This  Is  shown  In  Figure  3-8(D). 

Results  of  Image  reconstruction  experiments  applying  phase 
retrieval  to  the  estimate  of  the  object's  energy  spectrum  are  shown  In 
Figure  3-9.  Figure  3-9 (A)  shows  the  averaged  energy  spectrum  (with  the 
DC  term  removed)  of  the  pupil -plane  speckle  Intensity  for  the  wlder- 
segment  Golay-6  array  shown  In  Figure  3-6,  for  which  N  ■  10,240 
Independent  realizations  of  speckle  Intensity  were  averaged.  Figure 
3-9 (B)  Is  an  estimate  of  the  Fourier  modulus  of  the  object's  brightness 
distribution  computed  by  taking  the  square  root  of  the  Fourier 
transform  of  Figure  3-9(A).  Negative  numbers,  resulting  from  noise 
associated  with  the  finite-average  approximation  to  ah  ensemble 
average,  were  set  to  zero  prior  to  taking  the  square  root.  Figure 
3-9(C)  Is  the  Image  produced  by  applying  the  Iterative  transform  phase- 
retrieval  algorithm  [3. 6-3. 6]  to  the  Fourier  modulus  data  contained  In 
Figure  3-4(B).  The  procedure  for  accomplishing  phase  retrieval 
Involved  applying  several  cycles  of  the  hybrid  Input-output  algorithm 
(using  beta  ■  0.7)  and  the  error  reduction  algorithm  until  the 
algorithm  appeared  to  stagnate.  The  object-domain  constraints  used 
were  non-negativity  (since  an  unspeckled,  or  Incoherent,  Image  Is  being 
reconstructed)  and  a  loose  support  constraint,  a  rectangle  half  the 
size  of  the  smallest  rectangle  enclosing  the  average  energy  spectrum  of 
the  observed  speckle  pattern.  The  object  Is  guaranteed  to  fit  within 
this  support  constraint  [3.16], 


FIGURE  3-9.  IMAGE  RECOVERY  USING  THE  GOLAY-6  APERTURE,  N  ■  10,240. 
(A)  Average  energy  spectrum  of  the  measured  speckle  patterns  (with  dc- 
term  removed);  (B)  estimate  of  the  Fourier  modulus  of  the  object 
brightness  distribution;  (C)  image  reconstructed  from  (B)  using  the 
iterative  transform  (phase  retrieval)  algorithm;  (D)  Wiener-like  filter 
for  the  Golay-6  aperture;  (E)  filtered  Fourier  modulus  estimate;  (F) 
image  reconstructed  from  (E);  (G)  original  incoherent  object;  H 
filtered,  incoherent  object;  (I)  results  of  filtering  (C) . 


35 


Note  that  the  recovered  Image  shown  In  Fig.  3-9 (C)  Is  very  noisy 
compared  with  the  original  Incoherent  object,  shown  In  Fig.  3-9(G), 
although  a  general  resemblance  of  ^he  object  has  been  recovered.  Noise 
In  this  reconstructed  Image  Is  due  to  the  fact  that  a  finite  (albeit 
large)  number  of  speckle  realizations  were  used  to  estimate  the  Fourier 
modulus.  To  reduce  these  noise  effects,  we  multiplied  the  FouHer 
modulus  estimate  shown  In  Fig.  3-9(8)  by  the  Wiener-like  filter  of  Eq. 
(3-11).  For  these  simulations,  the  energy  spectrum  of  the  object  was 
taken  to  be  an  angular  (spin)  average  over  the  squared  Fourier  modulus 
of  the  true  object.  The  noise  spectrum  was  approximated  by  a  constant, 
whose  value  was  obtained  by  averaging  the  squared  Fourier  modulus 
estimate  over  those  higher  spatial  frequencies  where  the  signal -to- 
nolse  ratio  was  less  than  one.  Figure  3-9(D)  shows  the  resulting 
Wiener  filter  used  for  this  example.  Figure  3-9(E)  shows  the  product 
of  the  filter  3-4(D)  with  the  original  Fourier  modulus  estimate  3-9(B). 

Figure  3-9(F)  shows  the  Image  reconstructed  from  the  filtered 
Fourier  modulus  estimate  3-9(E)  using  the  phase  retrieval  algorithm, 
Note  that  the  filter  has  significantly  Improved  the  quality  of  the 
reconstructed  Image  3-9(F)  over  that  In  3«9(C)  reconstructed  without 
filtering.  For  the  purposes  of  comparison,  the  original  object  3-9(G) 
was  passed  through  the  filter  3-9(0),  with  the  result  shown  In  3-9(H). 
The  Image  reconstructed  from  speckle  correlation  measurements,  shown  In 
Fig.  3-9(F),  compares  favorably  with  the  filtered  object  3«9(H), 
Indicating  good  performance  on  the  part  of  the  Iterative  transform 
algorithm.  Figure  3-9(1)  shows  the  result  of  applying  the  Wiener 
filter  to  the  reconstructed  Image  shown  In  Fig.  3-9(C).  It  appears,  at 
least  for  this  example,  that  filtering  followed  by  Image  reconstruction 
Is  somewhat  superior  to  Image  reconstruction  followed  by  filtering.  We 
might  expect  to  get  even  better  results  by  using  an  Improved  Wiener 
filter,  for  example,  by  using  a  better  estimate  of  the  object  power 
spectrum  or  by  using  Eqs.  (3-13)  and  (3-14). 
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One  way  to  evaluate  the  MTF-boostlng  properties  of  the  filter  of 
Eq.  (3-11)  Is  by  Inspection  of  the  filter,  which  Is  shown  In  Figure 
3-9 ( D) .  Notice  that  It  has  a  local  minimum  In  the  center  (at  zero 
spatial  frequency)  and  a  ring  of  local  maxima  at  a  higher  spatial 
frequency.  This  compensates,  In  part,  for  the  rapid  drop-off  of  the 
OTF  that  can  be  seen  In  Figure  3-7,  The  ratio  of  the  peak  value  of  the 
filter  to  the  zero-frequency  value  Is  3.38,  a  sizable  boosting  of  the 
OTF  at  that  spatial  frequency.  This  falls  short  of  a  complete 
compensation  due  to  the  noise  energy  spectrum  term  In  Eq.  (3-11).  For 
the  same  reason,  the  filter  drops  off  for  the  highest  spatial 
frequencies,  where  the  noise  dominates  the  signal. 

Another  way  to  evaluate  the  MTF-boostlng  properties  of  the  filter 
of  Eq.  (3-11)  Is  to  compare  the  Imaging  results  shown  In  Figure  3-9 
with  those  obtained  with  a  filled  collecting  aperture;  Figure  3-2  In 
Section  3.1  shows  the  results  of  Image  recovery  from  simulations  of 
Imaging  correlography  obtained  with  a  full  aperture,  where  the 
simulated  speckle  Intensity  data  were  filtered  by  a  square  aperture 
comprising  64  x  64  "detector"  pixels  fully  encompassing  the  sparse 
Golay  aperture  used  above.  (The  width  of  the  Golay  array  Is  only  55 
pixels.)  Except  for  the  form  of  filtering  used  to  mask  the  speckle 
measurement  data,  the  digital,  processing  steps  used  to  produce  each 
frame  of  Figure  3-2  Is  Identical  to  that  of  the  corresponding  frame  of 
Figure  3-9.  The  top  row  of  frames  pf  Figure  3-2  correspond  to  Image 
retrieval  with  a  full  aperture,  but  without  Wiener  filtering.  Note 
that  the  resulting  Image  3-2(C)  Is  noisy,  but  Is  significantly  better 
than  Its  sparse  array  counterpart  3-9(C).  The  filter  shown  In  3-2 (D) 
is  that  prescribed  by  Eq.  (3-11)  with  the  OTF  given  by  the  auto¬ 
correlation  of  the  filled,  square  aperture.  Figure  3-2(F)  shows  the 
Image  recovered  from  the  Wiener-filtered  Fourier  modulus  3-2(E)  for  the 
filled  aperture.  Comparing  Figures  3-2(F)  and  3-9(F)  Indicates  that 
most  of  the  key  features  of  the  object  recovered  in  the  fllled-aperture 
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case  were  also  recovered  with  the  sparse,  Golay-6  aperture  case. 
However,  some  of  the  finer  details  of  the  object  recovered  In  the.  full 
aperture  case  were  smoothed  over  In  the  Golay  aperture  reconstruction. 
This  loss  of  resolution  for  the  sparse-aperture  case  Is  the  result  of  a 
smaller  OTF(Au)  value  (1.6.,  a  lower  redundancy),  and  hence  a  lower 
slgnal-to-nolse  ratio,  for  larger  spatial  frequencies. 

The  results  of  this  section  demonstrate  the  possibility  of 
recovering  Images  from  non Imaged  (far-fleld)  laser  speckle  patterns 
observed  with  sparse  arrays  of  Intensity  detectors.  The  Images 
obtained  using  a  combination  of  a  Wiener-filtered  speckle 
autocovariance  together  with  the  Iterative  transform  phase  retrieval 
algorithm  show  marked  Improvement  over  those  obtained  without 
filtering.  The  fact  that  the  Image  In  Fig.  3-9(F),  constructed  with 
sparse  arrays  of  detectors,  approaches  the  quality  of  the  full -aperture 
Image  shown  In  Figure  3-2(F)  suggests  that  the  MTF  boosting  filtering 
Is  successful  In  removing  Image  artifacts  due  to  the  sparse  collecting 
aperture. 

1  ■  •  ■  i 

Figure  3-10  shows  similar  results  for  the  sparse  aperture,  but  with 
only  N  -  1024  realizations  averaged.  As  In  the  fllled-aperture  case, 
the  loss  In  resolution  due  to  a  lower  signal -to-nolse  ratio  Is  evident. 

3.3  SIGNAL-T0-N0ISE  RATIO  AND  RESOLUTION 

Up  to  this  point,  we  have  alluded  to  the  fact  that  the  slgnal-to- 
nolse  ratio  and  the  quality  of  the  reconstructed  Images  In  Imaging 
correlography  Increases  with  the  number  N  of  Independently  observed 
speckle  patterns.  More  to  the  point,  the  error  In  the  speckle 
autocovariance,  and  so  tne  Fourier  modulus  estimate,  will  Improve  as 
the  number  of  speckle  measurements  Increase*,  whether  these  speckle 
measurements  arise  from  additional  speckle  pattern  realizations 
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FIGURE  3-10.  IMAGE  RECOVERY  USING  GOLAY-6  APERTURE,  N  -  1024.  (A) 
Average  eneray  spectrum  of  the  measured  speckle  patterns  (with  dc  term 
removed);  (8)  estimate  of  the  Fourier  modulus  of  the  object  brightness 
distribution;  (C)  Wiener-like  filter  for  the  Golay-6  aperture;  ( D) 
filtered  Fourier  modulus  estimate.*  (E)  Image  reconstructed  from  (D); 
(F)  original  incoherent  object;  (G)  filtered,  Incoherent  object. 


(snapshots)  or  from  an  increased  redundancy  in  the  OTF  of  the 
collecting  aperture.  This  flexibility,  in  choosing  between  number  of 
snapshots  N  and  collecting  array  redundancy,  can  be  better  appreciated 
by  considering  the  signal -to-noise  ratio  (SNR)  of  the  autocovariance 
estimate  achieved  in  imaging  correlography.  Assuming  that  time- 
sequential  measurements  of  the  speckle  patterns  are  statistically 
Independent,  we  can  show  that  the  SNR  of  the  estimate  of  the  object's 
energy  spectrum  at  spatial  frequency  Au  provided  by  the  estimator  of 
Eq.  (3-7)  is  given  by  [3.17,  3.18] 


SNRc(Au;  N)  - UlftHU! — __ 

C  [3  +  14l/»(Au)lZ  +  3l^(Au)r]1/Z 


(3-15) 


where  N  is  the  number  of  independent  speckle  patterns  (snapshots) 
observed,  /»( Au)  ■  r(Au)/r(0)  is  the  complex  coherence  factor  for  the 
measured  speckle  field,  and 


L  -  L (Au)  •  N5  OTF (Au)  (3-16) 


Is  the  number  of  redundant  pairs  of  speckle  intensity  in  the  collecting 
aperture  measured  at  pixel  separation  Au.  In  the  above,  Ns  is  the 
number  of  Independent  samples  of  Intensity  (or  number  of  speckles) 
contained  in  the  measurement  aperture  H(u).  For  the  case  that  the 
noise  in  the  Fourier  modulus  estimate  is  dominated  by  statistical 
fluctuations  In  the  autocovariance  estimate  Itself  (not  by  photon  shot 
noise,  etc.),  Eq.  (3-15)  specifies  the  trade-off  between  array 
redundancy  L  and  number  of  speckle  snapshots  N  needed  to  keep  the  SNR 
of  the  estimate  at  an  acceptably  high  level.  Keeping  the  SNR  of  the 
speckle  autocovariance,  and  so  the  SNR  of  the  estimate  of  the  object's 
Fourier  modulus,  at  a  high  level  will  preserve  an  acceptable  quality  in 
the  image  recovered  using  the  iterative  transform  phase-retrieval 
algorithm. 
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Notice  from  Eq.  (3-15)  that  the  SNR  Increases  with 

(1)  Increasing  aperture  size  [for  a  filled  aperture  SNR  a  (L)1^2  a 
(Ns)^2  a  (aperture  area)1^2. 

(2)  Increasing  object  diameter  [SNR  a  (Ns)^2  a  (object  area)1^2] 
and 

(3)  Increasing  number  of  snapshots  [SNR  a  (N)^2]. 

For  high  spatial  frequencies!  l^(Au)r  «  1  and 

SNRC(AUJ  N)  «  mpsrmr  \ftM\2/l T  .  (3-17) 

These  higher  spatial  frequencies  are  of  Interest  In  obtaining 
resolvable  detail  In  the  reconstructed  image. 

In  order  to  test  the  accuracy  of  these  SNR  expressions,  we 
computer-simulated  speckle  frames,  performed  the  correlography 
averaging,  and  determined  the  error  In  the  estimate  of  l/il2.  Taking 
the  Fourier  transform  of  the  data  averaged  according  to  Eq.  (3-3) 
yields  an  OTF-welghted  estimate  of  the  power  spectrum  of  the  Incoherent 
object,  as  given  by  Eq.  (3-10).  Normalising  that  to  unity  at  zero 

spatial  frequency  yields  an  estimate  of  OTF(Au)  l/i(Au)l2.  The  variance 
the  estimate  of  l/il2,  from  Eq.  (3-15),  Is 

Var{l/u(Au)  I2}  -  [3  +  14l>*(Au)  I2  +  3l/a(Au)  l4]/[NNs  OTF(Au)]  .  (3-18) 
Therefore  the  variance  of  the  estimate  of  OTF(Au)  l/*(Au)l  Is 

Var{0TF(Au)  lJ(Au) l2>  ■  OTF(Au)  [3  +  14l/a(Au) I2  +  3l/i(Au)l4]/(NN5)  . 

(3-19) 
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In  order  to  compare  theory  with  simulation  results,  It  Is  necessary 
to  compute  the  statistics  of  the  simulated  data  over  large  areas  In  the 
frequency  domain.  However,  the  variability  of  OTF(Au)  over  such  an 
area  would  confuse  the  results.  Hence  we  also  looked  at  an  estimate  of 
rOTF(Au)  \p\2,  which  has  a  variance 

Var{(OTFTiuT"  l/i(Au)l2}  ■  OTF(Au)  Var{lji(Au) I2} 

■  [3  +  14 I/s (Au)  I2  +  3l/i(Au)  I 4]/(NNs)  (3-20) 

which  Is  Independent  of  OTF(Au).  Thus  we  considered  the  two  absolute 
errors 

\  '■ 

ew(Au)  »  OTF(Au)  lJ(Au)l2  -  OTF(Au)  l/(Au)l2  '(3-21) 

wi^h  the  natural  OTF  weighting,  with  variance  given  by  Eq.  (3-19),,  and 

*  ,  •  '  ’  1  ‘M 

eu(Au)  JWTSuT  lJ(Au)l2  -  ^6tf(Aur  l/s(Au) I2  (3-22) 

with  JCTT  weighting,  with  variance  given  by  Eq,  (3-20),  which  Is  not 
weighted  by  the  OTF.  This  flJTT -weighted  data  was  obtained  by  dividing 
the  natural ly-OTF-welghted  data  by  fOTOTTT  [the  result  was  set  to 
zero  where  OTF (Au) »0] • 

Figure  3-11  shows  these  two  absolute  errors,  for  the  fllled- 
aperture  case  described  In  Section  3,1,  The  middle-grey  areas  have 
zero  error,  the  lighter  areas  have  positive  error,  and  the  darker  areas 
have  negative  error.  Figures  3-12  and  3-13  show  the  same  thing  for  the 
Golay  aperture  ckse  described  In  Section  3.2.  Note  that  for  the  OTF- 
welghted  case,  the  error  Is  maximum  near  Au»0  (at  the  center)  and  falls 
to  zero  at  the  edges  where  OTF(Au)  «  0,  as  predicted  by  Eq.  (3-19). 
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FIGURE  3-12.  SQUARED  FOURIER  MODULUS  ERROR  FOR  GOLAY-6  APERTURE 
(N  »  128). 


44 


FIGURE  3-13.  SQUARED  FOURIER  MODULUS  ERROR  FOR  GOLAY-6  APERTURE 
(N  •  1024,  10,240). 


For  the  f5TP -weighted  ("unweighted")  case  the  errors  are  more 
uniformly  distributed  across  the  frequency  plane,  as  predicted  by  Eq. 
(3-20).  Because  we  normalize  the  estimated  data  to  unity  at  Au-0,  by 
definition  we  made  the  absolute  error  equal  to  zero  at  Au>*0.  It  Is 
interesting  to  note  that  the  correlation  distance  of  the  absolute  error 
appears  to  decrease  as  the  number,  N  (K  In  the  figures),  of, frames 
averaged,  Increases.  This  effect  Is  presently  not  understood. 

The  averaging  to  compute  the  statistics  of  the  error  was  done  over 
a  32  x  32-plxels  square  area  shown  In  Figure  3-13 (A)  for  the  Gol ay- 
aperture  case  and  In  the  same-sized  square  In  the  corner  of  the  square 
fllled-aperture  case.  In  these  areas  J6Tf(au)  “  0.20  for  the  filled 
aperture  and  ~  0.10  for  the  Golay  aperture.  In  both  cases,  since  the 
object  fits  within  a  rectangle  of  size  40  x  60  pixels  embedded  in  a  12B 
x  128  array,  the  number  of  samples  per  speckle  In  the  aperture  plane  Is 
(128  x  128)/ (40  x  60)  -  6.83.  The  areas  of  the  filled  and  Golay 
apertures  were  measured  to  be  64  x  64  -  4096  and  822  pixels, 
respectively.  Consequently,  the  value  of  Ns,  the  number  of  speckles  In 
the  aperture  Is  4096/6.83  ■  600  for  the  filled  aperture  and  822/6.83  ■ 
120  for  the  Golay  aperture.  These  areas  of  Integration  were  chosen  to 
be  at  large  lAul  for  which  l/t(Au)l  «  1,  so  that  the  theoretical 
variance  expressions  simplify  to 


Var[0TF(Au)  l?(Au) lZ]  *  3  0TF(Au)/(NNs)  (3-23) 

and 

VarCJTJTFTSuT  l?(Au)l2]  *  3/(NNs)  .  (3-24) 

In  addition,  for  the  JOTF'  l/il2  case,  averages  of  the  statistics 
were  taken  over  the  entire  array  In  order  to  obtain  better  statistics. 
In  this  case  It  was  also  assumed  that  l/il  «  1  enabling  us  to  use  Eqs. 
(3-23)  and  (3-24).  Although  It  Is  not  true  for  small  values  of  lAul, 
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these  areas  are  relatively  small  compared  with  the  total  array  size,* 
consequently,  when  averaging  over  the  entire  array  the  assumption  that 
1^1  «  1  Is  reasonable. 

Table  3-1  compares  the  theoretical  expressions  for  the  noise 
variance  given  by  Eqs.  (3-23)  and  (3-24)  with  the  measured  variance  of 
the  simulated  data.  Individual  values  do  not  agree  very  well  because 
It  was  not  possible  to  Integrate  over  large  enough  areas  to  get  good 
statistics.  The  last  column,  the  ratio  of  the  measured  variance  to  the 
theoretical  variance,  shows  that  for  roughly  half  the  cases  the 
measured  variance  exceeded  the  theoretical  variance,  and  for  the  other 
half  It  was  less.  Thus  on  average  the  measured  noise  variance  roughly 
agrees  with  the  theoretical  expression,  giving  confirmation  of  the 
theory. 

Next  we  compute  the  number  of  frames  required  to  achieve  a  given 
resolution  for  a  particular  Imaging  scenario.  Figure  3-14  [3.19]  shows 
the  spin  (angularly)  averaged  l/i(Au)l2  for  the  satellite  object  at  a 
finer  resolution  than  for  the  version  of  the  linage  shown  earlier  In 
this  Section.  For  this  version  the  object  fit  within  a  128  x  128  array 
embedded  In  a  256  x  256  array,  and  Its  Fourier  transform  was  not 
weighted  by  an  OTF  function.  Hence  Its  resolution  was  about  4  times 
better  In  each  dimension  than  the  Image  shown  In  Figure  3-2(G).  The 
value  of  I Aul  at  the  highest  spatial  frequency  shown  In  Figure  3-14  Is 
therefore  equivalent  to  collecting  an  array  of  about  128  x  128 
speckles.  Thus  an  aperture  of  size  32  x  32  speckles  would  achieve  a 
resolution  equivalent  to  Au  ■  0.25  on  this  chart,  at  which  |/*l2  » 
0.012.  At  half  this  spatial  frequency  (l.e.,  Au  ■  0.125  on  this  plot) 
one  would  have  l^il2  ~  0.02.  Achieving  resolution  at  one-half  the 
highest  spatial  frequency  passed  by  the  aperture  Is  a  reasonable  goal, 
since  for  that  spatial  frequency  OTF(Au)  -  0.5,  whereas  for  higher 
spatial  frequencies  the  SNR  given  by  Eq.  (3-17)  drops  off  rapidly 
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Table  3-1 

, 

Comparison  of  Noise  Variance  Theory 
with  Simulations 


For  OTF(Au)  l/»(Au)l2,  32  x  32  regions  where  l/»l2  «  li 


Aperture' 

N 

'•  : 2  Theory  .  p 

i  ,  - .  ■ 

2 

Measured 

Ratio 

Filled 

1,024 

.10  E-5 

.40  E-6 

■  • 0.4 

Filled 

1,024 

.10  E-5 

.62  E-6 

0.6  , 

Filled 

10,000 

.10  E-6 

.19  E-6 

1.9 

Filled 

10,000 

.10  E-6 

.25  E-6 

2.5 

■  i  i 

Golay 

128 

.20  E-4 

.37  E-4 

1.8 

Golay  . 

1,024 

.25  E-5 

.13  E-5 

0.5 

Golay 

10,240 

.25  E-6 

.12  E-6 

0.5 

For  JOTF(Au)  ‘ 

■li(4u)l2, 

ehti re  array,  assume  t/»! 

2  «  1 1 

■  •.  r  • 

Aperture  _Jj _  g2  Theory  g2  Measured  1  Ratio 


Filled 

1,024 

.3  E-5 

.26.  E-5- 

0.9 

Filled 

10,000 

.3  E-6 

.  .76  E-6  . 

2,5 

Golay 

128 

.20  E-3 

.26  E-3 

1.3 

Golay 

1,024 

.24  E-4 

,15  E-4 

0.6 

Golay 

■  10,240 

.24  E-5 

,13  E-5  ' 

0.5 
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2 

because  both  OTF(Au)  and  l^(Au)l  are  decreasing  rapidly.  Plugging 

these  numbers  (N.  ■  32  x  32,  OTF  ■  0.5,  till2  ■  0.02)  Into  Eq.  (3-17) 

5  1 

yields  SNRC  *  0.26N  '  .  Therefore  In  order  to  achieve  a  given  SNRC  for 
this  example,  the  required  number  of  Independent  frames  Is  N  ?  15  SNR2. 
If,  for  example,  a  faithful  Image  requires  SNRC  -  10,  then  1,500  frames 
would  be  required  to  achieve  a  faithful  Image.  Recall  that,  for  this 
example,  since  32  x  32  speckles  were  assumed  to  be  In  the  filled 
aperture,  and  1/2  the  maximum  spatial  frequency  was  reconstructed,  then 
this  Image  would  have  a  space-bandwidth  product  of  16  x  16  resolution 
elements  (e.g.:  20  cm  resolution  for  an  object  of  width  3.2m  In 
diameter). 

If  the  array  size  were  doubled  In  each  dimension  In  order  to 
achieve  10  cm  resolution,  then  Ns  ■  64  x  64,  keep  OTF (Au)  -  0.5,  from 
Figure  3-14  !/»l 2  ~  0.012,  and  SNRC  -  0.31  N1^2;  and  one  needs  N  ~  10 
SNRC.  For  $NRC  ■  10,  this  requires  N  -  1,000  frames.  Thus  from  this 
example,  we  see  that  as  the  array  size  and  resolution  Increase,  the 
number  of  frames  can  decrease.  However,  this  Is  not  always  the  case: 
It  depends  upon  whether  the  product  of  the  aperture  length  with 
I ^ (Au) 1 2  Increases  or  decreases  as  the  aperture  length  and  Au  Increase 
together;  this  depends  upon  the  characteristics  of  the  target.  An 
object  consisting  of  a  small  collection  of  point-like  scatterers  or 
dominated  by  glints  will  have  a  l^(Au)l2  that  falls  off  much  more 
slowly  than  for  a  uniform,  very  extended  object.  Consequently  the 
point- like  objects  require  far  fewer  frames  to  achieve  a  given 
resolution.  Note  also  that  for  the  example  of  the  larger  aperture,  In 
order  to  get  the  20  cm  resolution  one  could  be  required  to  gather  only 
N  «  4  SNR2,  or  400  frames  for  SNRC  ■  10,  versus  1,500  frames  required 
for  the  smaller  aperture.  Thus,  we  see  the  trade-off  between  the  array 
area  (In  particular  Its  redundancy)  and  the  number  of  frames  In  order 
to  achieve  a  given  SNRC  at  a  given  spatial  frequency:  It  Is  the 
product  NN$  OTF(Au)  that  matters. 
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FIGURE  3-14.  SQUARED  FOURIER  MODULUS  FOR  P72-2  SATELLITE  MODEL  (FROM 
REF.  3.19).  ..  .  .  ..... 


SO 


The  above  noise  analysis  Included  only  the  effects  of  approximating 
the  ensemble  average  by  averaging  over  a  finite  number  of  realizations. 
The  variance  of  the  noise  that  Includes  both  finite  averaging  noise  and 
photon  noise  Is  given  by  [NNS  OTF(Au)]"1  times  [3.18] 


where  the  first  set  of  terms  Is  due  to  finite  averaging  and  the  second 
Is  due  to  photon  noise,  M  Is  the  number  of  detectors  (pixels)  per 
speckle,  and  <n>  Is  the  mean  number  of  photons  per  detector.  Thus  for 
l/il2  «  1,  the  finite  averaging  noise  variance  Is  proportional  to  3, 
while  the  photon  noise  variance  Is  proportional  to  4/(M  <n>)  +  1/ (M 
<n>2).  For  M  -  4,  the  two  noise  variances  are  equal  for  <n>  ■  1/2 
photon  per  detector  or  M  <n>  ■  2  photons  per  speckle.  Consequently, 
Independent  of  the  array  redundancy  and  the  number  of  realizations, 
photon  noise  will  be  negligible  as  long  as  the  number  M  <n>  of  photons 
per  speckle  Is  much  greater  than  two.  '• 

In  summary,  we  have  demonstrated  via  computer  simulations  that  It 
Is  possible  In  principle  to  recover  an  Incoherent  Image  of  a  laser- 
illuminated  object  from  multiple  realizations  of  detected  speckle 
Intensities  collected  over  sparse  arrays.  This  would  permit  the 
reconstruction  of  fine-resolution  Images  despite  phase  errors  due  to 
atmospheric  turbulence.  The  expressions  for  signal -to-nolse  ratio  as  a 
function  of  spatial  frequency,  array  redundancy,  and  number  of  speckle 
realizations  shows  that  large  amounts  of  array  redundancy  and/or  large 
numbers  of  speckle  realizations  are  required  to  reconstruct  an  Image  of 
large  space-bandwidth  product. 

Sections  3.2  and  3.3  are  an  expansion  of  References  3.20  to  3,22. 
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3.4  WIENER  FILTERING 


‘  *  .  |  -*»  ■  ■  ' 1 

The  Wlener-Hel strom  filters  given  by  Eqs.  (3-11)  and  (3-14)  y*ere 
shown  to  Improve  Image  quality;  however,  they  are  probably  not  the 
optimum  linear  filters.  In  this  section,  we  describe  two  different 
filters  that  could  yield  better  results  and  a  different  formulation  of 
the  problem  that  would  lead  to  a  different  filter.  Since  none  of  these 
were  Implemented  and  proven,  they  should  be  considered  speculative  at 
this  point. 

3.4.1  Recursive  Wiener  Filter 

In  the  Wlener-Hel strom  filters  given  by  Eqs.  (3-11)  and  (3-14).  we 
must  know  both  an  average  power  spectrum  of  the  signal  (estimated  by 
in2  and  IT!4,  respectively)  and  the  power  spectrum  of  the  noise, 
En(Au).  The  noise  can  often  be  determined  by  measuring  the  signal  plus 
noise  In  a  region  of  spatial -frequency  space  where  the  signal  Is  small. 
In  practice,  estimation  of  the  power  spectrum  of  the  object  Is  a  bigger 
problem.  Use  of  spin-averages  of  an  ensemble  of  Images  from  the  class 
of  objects  of  Interest  Is  one  approach.  Another  approach  Is  to  use  the 
realization  of  the  given  data.  This  can  be  done  by  the  method  that 
follows. 

For  simplicity,  first  consider  the  conventional  Imaging  model: 

,  ■  M 

t 

0  -  *  *  90  + -n  (3-25) 

where  g0  ■  the  Ideal  Image 

s  ■  the  psf 

n  ■  noise 

g  -  the  measured  Image 
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which,  In  the  Fourier  domain  Is 


G  ■  SfiQ  +  N  .  (3-26) 

Assuming  g0  and  n  are  Independent  stochastic  processes,  both  zero-mean 
Gausslans,  then  the  1e«st  mean-squared  error  linear  estimator  for  gQ  is 

g0  *  7-1  [gJ  where  [3a4] 

(3-27) 


(3-28) 

where  INI4  and  IGQI  are  the  power  spectral  densities  of  the  object  and 
the  noise,  respectively.  Helptrom  notes  that  although  the  Images 
generally  will  not  satisfy  the  statistical  assumption  stated  above, 
"Nevertheless...  It  can  be  expected  to  be  effective  when  applied  to 
these  Images  as  well."  So  the  Wlener-Hel strom  filter  Is  not  optimum, 
but  It  should  be  effective  and  It  Is  very  simple  to  use. 

A  big  problem  with  the  Wlener-Hel strom  filter  Is  that,  although  S 
Is  often  known  very  well  from  the  system  parameters  and  the  expected 
INI2  can  be  analyzed  ahead  of  time,  the  object  power  spectrum,  IGQI2  Is 
usually  unknown.  Helstrom  suggested  using  a  constant  for  IGfll2. 
However,  this  Is  a  very  bad  assumption  since  It  Is  well  known  that  the 
power  spectra  of  real,  nonnegative  objects  has  a  large  peak  at  low 
spatial  frequencies  and  then  drops  off  at  higher  frequencies,  typically 
at  a  rate  proportional  to  (spatial  frequency)*1. 
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t  , 

A  new  method  Is  to  use  the  data  Itself  as  the  estimate  of  the 
object  power  spectrum.  Namely,  use 

IG0I2  -  1 6 ! 2  .  (3-29) 

However,  as  seen  from  Eq.  (3-26),  this  may  be  a  poor  estimate.  Using 
this  estimate  In  Eqs.  (3-27)  and  (3-28)  yields  the  Wiener-filtered 
Fourier  transform 


where 


S* 

ISI2  +  I  Nr/ 1 61* 


(3-30) 

(3-31) 


Now  IGjI2  Is  a  better  estimate  of  IG0I2  than  1GI2  Is.  Therefore  It 
makes  sense  to  form  a  new  Wiener  filter 


S* 


*  +  INI* 


ISI 


INI  /IGjI' 


(3-32) 


and  apply  this  to  the  data.  This  supposedly  will  yield  a  better 
estimate  still.  Thus  we  can  get  successively  better  approximations  of 
the  Wiener  filter,  Wn,  and  of  the  object's  power  spectrum,  IGn+1l2  “ 
IWnGI2,  via 


W. 


S» 

ISI*  +  INI2/IGnlz 


(3-33) 


The  recursive  estimation  of  the  object's  power  spectrum  will 
converge  to  a  stable  solution  when 
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2 


(3-34) 


IG.I2  ■  I Gn+1 1 2  ■  IG/ 


—SIS, _ 

+  IN!Z/IGnl 


1 


Solving  for  1 6W 1 2 1 

IGJ2  -  0 
or 

IGJ4  ISI4  +  IGJ2  ISI2  (2IMI2  +  IGI2)  +  INI4  -  0  (3-35) 

which  simplifies  to 

IGJ2  "  (2ISI2)"1  (lGI2  -  2INI 2  *  HgI4  -  4INI5  IGI2  ]  .  (3-36) 

Note  that  for  IGI2  »  INI2,  IGJ2  »pprdache*  1 6 1 2/ 1 S 1 2  for  the 
positive-square-root  solution,  the  expected  result.  In  this  same 
limit,  the  negative  square  root  yields  zercv  and  so  only  the  positive 
square  root  should  be  used  in  Eq.  (3-36).  Plugging  this  Into  Eq. 
(3-33)  yields 


S*[lGI2  -  2INI2  +  JlGI4  -  41 Nl 2  IGI2  ] 
ISI2  (lGI2  +  JlGI4  -  4 1 N I ^  IGI2  ) 


(3-37a) 


S+[l  +  Jl  -  4INI2/IGI2  ' 
2 1 S I  ^ 


(3-37b) 


Note  that  for  INI2/IGI2  «  1, 

W.  *  [l  -  -^4]  (3-38) 

isr-l  IGIZJ 

and  as  INI2/IGI2  4  (i/4)“, 

W,,  ♦  1  •  (3-39) 

2  isr 

For  INI2/IGI2  >  1/4,  Eq.  (3-37)  Is  invalid.  By  performing  the  recursive 
calculation  of  Eq.  (3-33)  for  this  situation,  It  was  found  that  IGn I2 
became  progressively  smaller,  so  that  for  INI2/! G I2  >  1/4, 

IGJ2  0  and  Wfl  •»  o  .  (3-40) 

It  has  not  yet  beon  assessed  which  of  the  Wiener  filters  discussed 
above  Is  best. 

Conversion  of  these  results  to  the  case  of  Imaging  correlography  Is 
accomplished  by  replacing  G  In  the  equations  above  by  Cj(Au;N). 
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3.4.2  Iterative  Nonlinear  Filter 


A  completely  different  approach  to  filtering  is  as  follows.  Suppose 
that  gQ(x)  and  g (x)  of  Eq.  (3-25)  are  spatially  limited  to  some  support 
defined  by 


x«fl0(x)  >  0 
xtgQ(x)  ■  0 


(3-41) 


for  a  nonnegative  g0(x)  and  g(x).  Further  it  Is  assumed  that  the 
support,  $.(u),  of  S(u)t 

5 


1  ,  u:$(u)  >  0 
,0  ,  u:S(u)  ■  0 


(3-42) 


is  known.  Then  we  can  use  the  Iterative  transform  algorithm  to  better 
estimate  gQ(x)*s(x)  and  G0(u)  S(u)  by  iteratively  setting  successive 
estimate',  of  G0(u)  S(u)  to  zero  wherever  $s(u)  -  0  and  setting 
successive  estimates  of  gQ(x)*s(x)  to  zero  whenever  it  is  negative  or 
where  g$(x)  Is  zero.  Since  functions  satisfying  these  conditions  in 
both  domains  form  convex  sets,  this  error-reduction  algorithm  is  a 
projections-onto-convex-sets  (POCS)  algorithm  which,  by  Youla's  analysis 
[3.23],  has  strong  convergence  properties  (It  may  not  be  unique, 
though).  This  Iterative  would  reduce  the  noise  and  Is  an  alternative  to 
Wiener  filtering. 


3.4.3  Improved  Noise  Model 


The  new  methods  for  filtering  described  above  attempt  to  makeup  for 
the  lack  of  knowledge  of  the  power  spectrum  of  the  signal.  In  this 
section  we  point  out  that  additional  analysis  needs  to  be  performed  in 
order  to  properly  model  the  noise. 
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The  two  signal-plus-noise  models  assumed  for  the  filtering  actually 
performed  were  [see  Eq.  (3-12)]  OTF(Au)  r(Au)  +  N  and  OTF(Au)  IT(Au) I2  + 
Nc.  After  the  image  reconstruction  experiments  were  oerformed,  the 
signal-to-noise  ratio,  proportional  to  given  by  Eqs. 
(3-15) - (3-17) ,  was  analyzed.  Since  the  signal  is  proportional  to  OTF, 
this  implies  that  the  noise  is  also  proportional  to  1 OTF  1 .  This  fact 
was  verified  by  the  analysis  and  digital  experimental  described  in 
Section  3.3.  Thus  in  order  to  derive  an  optimum  linear  filter  for 
imaging  correlography  it  will  be  necessary,  in  future  research,  to 
derive  a  new  Wiener-type  filter  based  on  a  model  In  which  the  additive 
noise  is  weighted  by  fCTF’. 
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4.0  IMAGING  AND  CORRELOGRAPHY  WITH  A  MIXED  OBJECT 
4.1  INTRODUCTION 

,1 

In  the  previous  section  it  was  assumed,  for  the  multiple 
realizations  employed  In  Imaging  corral ography,  that  the  optical  fields 
reflected  by  the  object  are  diffuse,  zero  mean,  and  uncorrelated. 
However  this  will  not  be  true  when  the  Object  has  one  or  more  glint 
components.  In  this  section  we  analyze  the  case  In  which  the  object's 
reflectivity  contains  both  a  diffuse  component  and  a  glint,  or 
deterministic,  component,  Here  we  consider  not  only  Imaging 
correl ography,  but  other  Imaging  modalities  as  well. 

j  .  i 

In  most  instances  we  model  the  Imaging  process  as  either  coherent 
or  Incoherent.  Froiii  a  coherent  system  one  can  get  an  Incoherent  Image 
either  by  (1)  noncoherently  averaging  the  Intensities  of  many  coherent 
Images  or  (2)  by  heterodyne  Interferometry  with  multiple  realizations 
or  (3)  by  Imaging  correl ography,  which  Involves  first  averaging  the 
autocorrelations  of  the  Intensities  of  the  aperture-plane  fields  to 
estimate  the  energy  spectrum  of  the  Incoherent  object,  then 
reconstructing  an  Image  by  phase  retrieval.  In  the  analysis  of  the 
formation  of  an  Incoherent  Image  from  multiple  realizations  of  coherent 
data,  one  usually  assumes  that  the  object  Is  diffuse.  In  what  follows 
we  analyze  the  case  for  a  mixed  object,  l.e.  one  that  contains  both  a 
deterministic  component  and  a  diffuse,  variable  component.  All  three 
modes  of  obtaining  an  Incoherent  Image  from  coherent  data  of  a  mixed 
object  will  be  analysed. 
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4.2  MIXED  OBJECT  MODEL 

We  assume  that  the  object  Is  coherently  Illuminated  and  that  Hi 
complex  amplitude  reflectivity  consists  of  two  components,  a 
deterministic  component  g(x)  (where  x  is  a  2-D  spatial  or  angular 
coordinate)  and  a  diffuse  component,  dn(x)Y 

f„(x)  -  g(x)  +  d^(x)  .  (4-1) 

The  subscript  n  Indicates  the  realization  number,  where  dn(x)  Is 
assumed  to  have  different  realizations  of  phase,  due  to  a  rough 
surface-height  distribution,  as  the  look  angle  changes  slightly.  The 
underlying  Incoherent  object,  l.e.  the  object  Intensity  reflectivity 
had  It  been  Illuminated  by  Incoherent  light,  of  the  diffuse  component 
Is  given  by 

dj(x)  -<ldn(x)l2>n  (4-2) 

where  <*>n  denotes  an  ensemble  average.  We  assume  that  dn(x)  Is  zero 
mean  and  spatially  uncorrelated: 

<dn(x)>n  -  0  (4-3a) 

and 

<dn(x1)  d*(x2)>n  -  dj(Xj)  5(xj  -  x2)  .  (4-3b) 

In  the  aperture  plane,  In  the  far-fleld  relative  to  the  object,  the 
complex  amplitude  of  the  backseat tered  radiation  Is  given  by  the 
Fourier  transform  of  the  coherent  object, 

F„(u)  -  G(u)  +  D„(u)  (4-4) 
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where  the  functions  represented  by  uppercase  letters  are  the  Fourier 
transforms  of  the  functions  represented  by  the  corresponding  lowercase 
letters,  and  u  Is  a  2-D  coordinate  In  the  aperture  plane.  The 
component  G(u)  Is  deterministic.  The  diffuse  component  has  the 
following  properties! 

<Dn(u)>n  -  0  (4-5) 

and 

<Dn(Uj)  Dj(u2)>n  -  rQ( ^  -  u2)  -  TD(Au)  (4-6) 

where  superscript  *  denotes  complex  conjugate,  Au  ■  Uj  -  u2,  and,  by 
the  van  Clttert-Zernlke  theorem, 

r0(4u)  ■  7tdj(x)]  («-7) 

where  7  denotes  Fourier  transformation.  I Dn (u) I  would  be  a  fully- 
developed  speckle  pattern  (with  negative  exponential  point  statistics) 
and  Dn(u)  Is  circular-complex  Gaussian  (ccg)  distributed. 

An  example  of  an  object  satisfying  the  assumptions  above  Is  one 
with  a  rough  component,  dn (x) ,  that  rotates  slightly  to  result  In 
different  realizations,  plus  an  unchanging  component  g(x).  The 
unchanging  component  could  be  a  part  that  does  not  rotate  or  It  could 
be  a  single  unresolved  corner  reflector  or  stable  glint. 

4.3  NONCOHERENT  AVERAGING  OF  COHERENT  IMAGES 

A  single  realization  of  a  coherent  Image  of  the  object  through  an 
aperture  A(u),  with  coherent  Impulse  response  a(x)  ■  T"*[A(u)],  Is 
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a(x)  *  fn(x)  -  a(x)  *  g(x)  +  a(x)  *  dn(x)  ,  (4-8) 

whsre  *  denotes  convolution.  The  average  image  Intensity  Is 

<l«(x)  *  fn>x)|2*n  “  *  0M|Z  +  ‘‘■•W  *  <liiW|2>n 

♦  [•*(«)*  9*<x)]  [«(x)  *  <d„(x)>„]  +  c.c.  (4-9) 

where  c.c.  denotes  complex  conjugate.  The  second  term  Is 

<l*(x)  *  dn(x)l2>n  »  <lj  atxj)  dn(x  -  Xj)  dxjl2^ 

■  J  I  a(xj)  a*(x2)  <dn(x  -  Xj)  dj(x  -  x2)>n  dxJ  dx2 

-  J  la(x1)l2  dj(x  -  xt) 

-  Ia(x)l2  *  dj(x)  “  •  (4-10) 

where  use  was  made  of  Eq.  (4-3b).  Using  this  and  Eq.  (4-3a)„  Eq.  (4-9) 
slmplles  to 

<l«(x)  *  fn(x)  l2^  -  la(x)  *  g(x)l2  +  la(x)l2  *  dj(x)  .  (4-11) 

'/  ; 

That  Is,  In  the  noncoherent^  averaged  Intensity  image  of  a  mixed 
object,  the  deterministic  component  Images  as  In  a  coherent  system 
(l.e.,  convolution  with  the  coherent  Impulse  response,  then  modulus 
squared)  while  the  diffuse  component  Images  as  In  an  Incoherent  systom 
(l.e.,  Incoherent  object  convolved  with  the  Incoherent  Impulse  response 
which  Is  the  squared  modulus  of  the  coherent  Impulse  response).  This 
result  holds  both  for  optical  Imaging  systems  and  for  microwave  SAR. 
More  about  this  Image  will  be  said  In  the  next  section. 
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4.4  HETERODYNE- INTERFEROMETRY  AVERAGING 

The  second  Imaging  mode  Is  to  measure  complex  fields  In  the 
aperture  plane  using  heterodyne  detection,  then  calculate  the  average 
coherence  function  for  multiple  Realizations.  This  differs  somewhat 
from  conventional  amplitude  Interferometry,  ,  since  (Ignoring  aperture 
effects  for  the  moment) 

rF<ul*  u2>  -  <Fn<ul) 

«  G(uj)  G*(u2)  +  <On(uj)  P*(u2)>n 

+  G*(uj)  <Dn(u2)>n  +  c,c* 

■  Gfuj)  G*(u2)  +  TptUj  -  u2)  ■  (4i-12) 

Is  not  stationary  due  to  the  presence  of  the  deterministic  term. 
Direct  Image  reconstruction  from  rp(Uj,  u2)  Is  not  possible  because  the 
two  terms  on  the  right-hand  side  of  Eq.  (4-12)  cannot  be  easily 
separated  out. 

Consider  the  use  of  spatial  Integration  or  averaging.  We  denote 

<  •  >s  “  J  *  (4-13) 

V 

We  also  write  <<*>n>s  «  <*>ns,  Including  the  aperture  functions 
explicitly  in  Eq.  (4-12)  ind  spatially  Integrating  the  ensemble 
average,  for  a  given  vdTue'of  Au  -  Uj  -  u2  (or  u2  -  Uj  -  Au),  yields 
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<a(uj)  a(u2)  rF(ux,  u2)>8 


<A(U|)  W  A(u2)  Fn(M2»ns 


-  «A(u1)  G(Uj)  A(u2)  G*(u2)>8  4  <A(Uj)  A(u2) rD(ur  -  u2)>5 
»  (AQ)f(AG) (Au)  +  S(Au)  rn(Au)  (4-14) 

'  .1.1*  X'  '  *  iltj-i 

where  (AG)l(AG)(Au)  denotes  the  autocorrelation  of  A(u)  G(u)  evaluated 
at  Au,  and  S(Au)  Is  proportional  to  the  optical  transfer  function  due 
to  the  aperture  A(u): 


S(Au)  ■  AIA(Au)  «  <A(Uj)  A(u2)>s 

-  A0  OTF(AU)  (4-15) 

.  i 

where  OTF(Au)  Is  the  optical  transfer  function  for  the  aperture  A(u) 
and 

A0  ■  /  A2(u)  du  ,  (4-16) 

i 

Inverse  Fourier  transforming  Eq.  (4-14)  as  a  function  of  Au,  one 
gets  the  Image 

la(x)  *  q(x) I2  4  la(x)j2  *  dj(x)  ■  <|a(x)  *  fn(x)  (4-17) 

which  Is  Identical  to  the  noncoherent^  averaged  Intensity  Image  In 
Eq.  (4-11),  We  refer  to  this  .Image  as  the  Incoherent,  image  of  the 
mixed  object.  It  Is  not  clear  whether  there  exists  an  Incoherent 
object  'that  would  give  rise  to  such  an  Image,  but  for  convenience  we 
define  the  Incoherent  mixed  object,  fj(x),  by 
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(4-18) 


I  a  (x)  1 2  *  f  j  (x)  ■  <!  a  (x)  *  f  n  (x)  1 2>n  . 

The  Fourier  transform  of  the  Incoherent  image  of  the  Incoherent  mixed 
object  is 


S(Au)  Fj(Au)  -  (AG)# (AG) (Au)  +  S(Au)  rQ(Au)  (4-19) 

where  Fj(Au)  Is  the  Fourier  transform  of  the  Incoherent  mixed  object. 

From  Eqs.  (4-1),  (4-2),  and  (4-3a),  the  noncoherently  averaged 

object  Is 

<lfn(x)l2>n  -  lg(x) I2  +<ldn(x)l2>n 

+  g*(x)  <dn(x)>n  +  c.c. 

••  \  • 

-  Ig(x)l2  +  dj(x)  (4-20) 

An  Incoherent  Image  of  the  noncoherently  averaged  object  would  be 

I  a (x) 1 2  *  <lfn(x)l2>n  -  la(x)l2  *  I g (x) 1 2  +  la(x)l2  *  dj(x)  .  (4-21) 

Note  that  this  generally  differs  from  the  Incoherent  Image  of  the  mixed 
object  In  Eqs.  (4-11)  and  (4-17),  except  when  g(x)  Is  a  single  delta- 
function,  In  which  case  they  are  the  same,  that  Is  la(x)l2  *  I g (x) 1 2  ■ 
I  a  (x)  *g(x)l2.  Consequently,  the  noncoherently  averaged  object  Is 
generally  different  from  the  Incoherent  mixed  object.  From  Eq.  (4-19), 
the  Incoherent  mixed  object  would  be  the  Inverse  Fourier  transform  of 
S”1 (Au)  C(AG)*(AG) (Au)]  +  Tq(Au).  However,  It  Is  not  clear  whether  the 
first  term  yields  a  reasonable  "object"  component  In  the  general  case. 
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4.5  CORRELOGRAPHY/ INTENSITY  INTERFEROMETRY 


In  conventional  correlography  (diffuse  component  only),  since  Dn(u) 
is  ccg,  the  Gaussian  moment  theorem  can  be  used  to  obtain  (Ignoring 
aperture  effects  for  the  moment) 


<IDn<ul>|2  IBn<“2>|2>„  -  l<0n(ul>  D«<“2>n|2  +  <,0n<ul> |2s”n  <IDn<u2> |2>n 

>  ir„(u1  -  u-,)!2  *T§  (4-22) 


where  T0  Is  the  average  Intensity,  so  that  the  desired  energy  spectrum 
of  dj(x)  can  be  obtained  from  the  measured  data: 

|rp<ui  -  u2)l2  -  <IDn(ul)l2  IDn(u2)12>n  "  T0  *  <4“23) 


Since  all  the  terms  In  Eq.  (4-23)  are  stationary,  spatial  averaging  can 
also  be  Included  in  Eq.  (4-23)  to  yield  the  samq  result. 

For  the  mixed-objeet  case,  represented  by  Eqs.  (4-1)  and  (4-4),  the 
expression  analogous  to  Eq.  (4-22)  has  sixteen  terms: 

<w*  wS 

,  »  IG(Uj) I2  |G(u2)I2  +  IG(Uj)I2  <IDn(u2)l2>n  +  IG(u2)!2  <IDn(ui)lZ>n 
+  IG(u1)I2[g*(u2)  <Dn(u2)>n  +  c.c.]  ■MG(u2)l2[G*(ni)  <®n(uI)>(|  +  c.c.] 

i  ■  .  .  .  •  H 

+  <IDn(u1)l2  IDn(u2)l2>n 
+  G*(u2)  <0n(u2)  I0n(u1)l2>n  +  c.c. 
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(4-24) 


+  G*(Uj)  <D„(u1)  IDn(u2) r2>n  +  c.c. 

+  G*(Uj)  G*(u2)  <Dn(u1)  0n(u2)>n  +  c.c. 

+  G*(u1)  G(u2)  <Dn(u.)  D*(u2)>n  +  c.c.  . 

Using  <On(Uj)>n  ■  0,  <Dn(u2)  IDn(U|)l^n  ■  0,  «nd  <Dn  <ul>  °n<u2>>n  ‘  »' 
this  simplifies  to 

<IFn(Ul)l2  IFn(u2)l2>n  -  IG(‘Uj)l2  IQ(u2)I2 

+  IG)^)!2  <IOn(u2)l2>„  +  IG(uz)l2  <IDn(Uj) 1 2>,| 

+  ir0(uj  -  u2)i2  ♦  <ion(uj)i2>  <iDn(u2)i2>„ 

♦  G*(uj)  G(Uj)  r0(Uj  -  u2)  +  c.c. 

*  <•  WSi  <IFn<u2>|Z>n  +  lrD<“l  ' 

♦  G*(Uj)  G(u2)  rD(Uj  -  u2)  +  c.c.  (4-25) 

where  we  have  used,  from  Eq.  (4-4), 

<|Fn(u)|2>n  -  <IG(u)  +  Dn(u)l2>n 

*  IG(u) I2  +  <|Dn(u) l2>n  +  <£*(u)  Dn(u)>n  c.c. 

«  IG(u)l2  +  <IDn(u) l2>n 

-  IG(u)l2  +  rD(0)  .  (4-26) 

Notice  that  Eq.  (4-25)  Is  not  stationary. 
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Now  consider  the  spatial  Integration  of  Eq.  (4-25) ,  in  which  we 
Include  explicit  aperture  functions} 

<A(uj)  IFn(u1) I2  A(u2)  IFn(u2)l2>ns 

•  <A(ut)  <IFn(Uj)l2>n  A(u2)  <IFn(u2)l2>B>, 

♦  <a(u1)  a(u2)  ir0<u1  -  u?) I *>* 

♦  <A(ut)  G*(Uj)  A(u2)  G(u2)  r0(Uj  -  U2)>4  +  C.C.  (4-27) 
-  (A  <IFnl2>n)«(A  <!Fnl2>n)(Au)  ♦  S(Au)  ir0(Au)l2 

+  [(AG)I(AG) (Au)]*  Tq(Au)  +  c.c.  (4-20) 

where  Au  «  u.  -  u2.  From  Eq.  (4-19)  we  have 

S2 (Au)  I F | (Au) 1 2  ■  I (AG)I(AG) (Au) i2  ♦  S2(Au)  irQ(Au) I2 

+  [(AG) I (AG) (Au)]*  S(Au)  rp(Au)  +  c.c.  (4-29) 

Inserting  this  Into  Eq-  (4-28)  yields 

<A(u,)  IFn(Ul)l2  A(u2) 

-  (A  <IFnlZ>n)«(A  <IFnl2>n)(Au} 

+  S (Au)  I F j (Au) 1 2  -  S" 1 (Au)  l(AG)*(AG)(Au|l2  (4-30) 
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where  we  define 


S’1 (Au) 


’1/S  (Au) 
.  0 


S (Au)  >  0 
S (Au)  -  0  . 


(4-31) 


As  In  conventional  correlography,  Eq.  (4-30)  relates  the  ensemble 
and  spatially  averaged  aperture-plane  Intensities  to  the  sum  of  an  OfF- 
welghted. power  spectrum  of  the  Incoherent  (mixed)  object  and  the  square 
of  the  average  Intensity  (or  more  precisely,  the  autocorrelation  of  the 
ensemble  averaged  Intensity).  However,  the  deterministic  Component 
does  not  add  to  the  bias  In  the  same  way  as  the  diffuse  component, 
requiring  the  subtraction  of  the  additional  (the  last)  term  In  Eq. 
(4-30). 

IGn(u)r  can  be  determined  as  follows.  For  u2  ■  Uj  »  u  (Au  -  0), 
Eq.  (4-25)  gives 

<IFn(u)l4>n  •  (<IFn(u)l*>n)^  +  rjj(0)  ♦  2IG(u)l2  rD(0)  (4-32) 


where  we  have  used  the  fact  that  rQ(0)  Is  real  valued.  Eqs.  (4-26)  and 
(4-32)  constitute  a  set  of  two  simultaneous  equations  In  two  unknowns, 
I G (u) 1 2  and  rD(0),  Their  solution  yields 

I G(u) 1 4  -  2(<IFn (u) l2>n)2  -  <IFn(u) l4>n  (4-33) 

and 

r0(Ci)  ■  <IFn(u)l2>n  -  IG(u)l2 

■  <IFn(u)l2>„  -  j2(<IFn(u)l2>n)2  -  <IFp(u) |V„  .  (4-34) 
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2  •».■■■ 

For  a  complete  solution  of  IFj(Au)l  from  Eq.  (4-30),  one  must 
determine  the  last  term  in  Eq.  (4-30)  from  the  data,  a  task  yet  to  be 
accomplished. 

Note  that  for  G(u)  •  0  (i.e.,  no  deterministic  component),  Eq. 
(4-30)  reduces  to  the  conventional  correlography  result,  and  for  Dn(u) 
■  0  (I.e.,  only  a  deterministic  component),  the  last  two  terms  (the 
interesting  ones)  In  Eq.  (4-30)  cancel,  leaving  a  useless  equality. 

Single-glint  case. 

Now  consider  the  specific  case  for  which  the  deterministic 
component  of  the  object  Is  a  single  delta  function  (a  glint  or  corner 
reflector): 


Then  we  have 


Q(x)  -  b  d(x  -  xQ)  . 


G(u)  -  b  exp(-12fruxQ)  , 


(4-35) 

(4-36) 


<lfn(x)l2>n 


fT(x) 


-  I bl z  d(x  -  xQ)  +  dj(x) 


(4-37) 


Fj(Au)  «  Ibr  exp(-12r  Au  .<p)  +  Tq(Au) 


(4-30) 


IFj(Au) !2  -  I b 1 4  +  ir0(Au)l2 


♦  Ibr  exp(12r  Au  xQ)  Fq(Au)  +  c.c.  ,  (4-39) 
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<IFn(u)l2>I1  -  Ibl2  ♦  <IOn(u)l2>n  •  Ibl2  +  r0(0) 

I  <IFnl2>„  (4-40) 

and 

<IF„(u)l4>n  •  Ibl4  ♦  4lbl2  rD(0)  +  zr|(0) 

■  <IF„i4>  (4-41) 

the  latter  two  equations  being  Independent  of  u.  For  this  object,  Eq. 
(4-30)  becomes 

^(“i)  IFn<ul>,Z  A<UZJ  IFntu2)|2>ns 

-  S(Au)  [(<IFnl2>n)2  +  IFj(Au)l2  -  Ibl4] 

■  S(Au)  [ I Fj(Au) 1 2  -  (<l Fnl 2>n)2  +  <IFnl\]  (4-42) 

where  we  used,  from  Eq.  (4-33), 

Ibl4  -  I6(u)l4  -  2(<IFnlZ>n)2  -  <!Fnl4>n  .  (4-43) 

2 

Therefore  we  can  solve  for  I Fj (Au) I  for  this  object  In  terms  of 
measurable  quantities; 

S  (Au)  I F  |  (Au)  I2  “  <A(llj)  IFn(Ul)l2  A(u2)  IF^Uj)!2^, 

.  S(AU)  [WF.I  V  -  <IFnl\] 

■  <A(uj)  IFn(Ul)l2  A(u2)  *Fn(u2)l2>ns  -  S (Au)  Varn(IFnl2)  , 

(4-44) 
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where  Varn( I I  )  is  the  variance  over  the  ensemble  of  intensity 
realizations.  The  incoherent  image  of  this  type  of  mixed  object  can  be 
reconstructed  from  Eq.  (4-44)  via  phase  retrieval  with  the  help  of  a 
nonnegativity  constraint. 

Note  that  this  differs  from  the  usual  Ijig  -  T2  averaging  done  for 
the  case  of  a  diffuse  object  only.  In  that  style  of  notation,  Eq. 
(4-44)  is  equivalent  to 


IFj(Au)!2  -  <1^2  -  7  +  .  (4-45) 

Since  for  a  fully-developed  speckle  pattern 


(4-46) 


this  estimator  Is  equivalent  to 

I Fj(Au)  I2  ■  cljlg  -  2T2  +T2>  »  <Ijl2  -  7^  (4-47) 

for  the  case  of  G(u)  ■  0.  That  is,  the  estimator  of  Eqs.  (4-44)  and 
(4-45)  gives  the  correct  result  whether  there  4  a  single  glint  or  no 
glint,  whereas  the  other  estimators,  such  as  Ijl2  -  T2  give  the  wrong 
answer  if  there  is  a  single  deterministic  glint.  The  presence  of  a 
glint  can  be  detected  by  computing  Eq.  (4-43). 

Notice  that  ensemble  averaging  alone  and  ensemble  averaging 
followed  by  spatial  averaging  yields  useful  information.  However,  we 
were  unable  to  find  useful  Information  in  spatial  averaging  alone  for 
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the  case  of  a  mixed  object.  (And  spatial  averaging  followed  by 
ensemble  averaging  Is  equivalent  to  ensemble  averaging  followed  by 
spatial  averagings  «#>n>s  ■  «*>s>n  for  the  quantities  analyzed.) 

This  type  of  Image  formation  would  be  unaffected  by  atmospheric 
turbulence  as  long  as  nonisoplanatlsm  and  scintillation  are  not  present 
since  only  aperture-plane  Intensities  are  measured.  However, 
atmospheric  turbulence  would  severely  limit  conventional  coherent 
imaging  or  heterodyne  Interferometry. 
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5.0  COHERENT  IMAGE  RECONSTRUCTION  FOR  OBJECTS  HAVING  GLINTS 
5.1  INTRODUCTION 

Ordinarily  when  we  image  a  space  object  with  a  large-aperture 
earth-bound  optical  imaging  sensor,  the  resolution  1*  very  poor  owing 
to  phase  aberrations  caused  by  the,  turbulent  atmosphere.  An  approach 
to  circumvent  this  problem,  Initially  studied  over  a  decade  ago,,  Is 
laser  correlography  [5.1].  It  Involves  the  Illumination  of  the  space 
object  by  a  coherent  laser  and  the  detection  of  the  backscattered 
intensity  pattern  In  the  aperture  plane.  This  speckled  Intensity 
pattern  Is  the  squared  modulus  of  the  Fourier  transform  of  the  complex¬ 
valued  object  reflectivity.  Ordinarily  this  allows  only  an 
autocorrelation  of  the  object,  not  an  image  of  the  object,  to  be 
computed. 

For  certain  favorable  object  geometries,  such  as  objects  having 
separated  parts,  an  Iterative  Fourier  transform  algorithm  has  been 
developed  that  retrieves  the  phase  of  the  Fourier  transform  of  the 
object  and  thereby  permits  a  diffraction-limited  comp! ex- valued  image 
to  be  reconstructed  [5.2].  Unfortunately  the  class  of  objects  for 
which  this  approach  currently  works  Is  too  limited.  An  approach  that 
works  for  general  objects  is  that  of  imaging  correlography  [5.3] 
described  in  Section  3.  By  averaging  over  the  autocorrelations  of  many 
aperture-plane  speckle  Intensity  patterns,  one  arrives  at  the  modulus 
of  the  Fourier  transform  of  the  incoherent  object  (the  object  had  it 
beon  Illuminated  by  an  incoherent  source  such  as  the  sun).  In  this 
case  the  ideal  Image  Is  real -valued  and  nonnegative,  and  an  Image  of  a 
general  object  can  be  reconstructed  from  the  Fourier  modulus  data  using 
the  iterative  transform  algorithm  employing  the  powerful  nonnegativity 
constraint.  This  concept  has  recently  been  verified  In  laboratory 
simulations  [5.4].  Unfortunately  In  many  circumstances  it  may  be 
impractical  to  collect  a  large  enough  number  of  snapshots  of  aperture- 
plane  speckle  patterns  to  arrive  at  a  good  statistical  estimate  of  the 
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Incoherent  Fourier  modulus.  For  this  reason  it  is  of  great  interest  to 
be  able  to  reconstruct  a  coherent  image  from  a  single  snapshot  of  data, 
despite  the  difficulty  mentioned  at  the  beginning  of  this  paragraph. 

♦ 

It  has  been  observed  that  satellites  frequently  have  strong  glint 
returns  (mirror-like  reflections  off  solar  panels,  for  example).  It 
has  long  been  known  that  if  there  is  a  Single  glint  sufficiently 
separated  from  the  Vest  of  the  object,  then  by  the  holographic  method 
[5.5]  one  can  easily  reconstruct  a  coherent  Image  from  a  single 
snapshot  of  data.  However,  such  Ideal  glints  would  be  relatively 
uncommon.  Glints  centered  on  the  object  or  multiple  glints  would 
prevent  the  use  of  the  holographic  approach.’ 

1  .  \  •  '  , 

In  this  section  we  describe  methods  developed  for  reconstructing  an 
Image  having  glints  from  a  single  realization  of  the  Intensity  of  the 
aperture-plane  speckle  pattern  from  a  coherently  Illuminated  object. 

In  Section  5.2  the  most  successful  method  we  developed  Is  described. 

It  consists  of  three  successlvel  algorithms.  In  Section  5.2  is  shown 
reconstructions  using  just  the  iterative  transform  algorithm.  A 
recursive  reconstruction  algorithm  that  gave  limited  success  Is 
described  in  Section  5.3.  In  Section  5.4  the  effect  of  a  large  glint 
on  the  quantization  error  of  the  measured  data  Is  analyzed. 

5.2  THREE-ALGOR I THM  METHOD 

t  . 

In  what  follows,  we  describe  an  approach  that  permits  a  high-  ♦ 

fidelity  Image  to  be  reconstructed  from  a  single  snapshot  of  aperture- 
plane  Intensity  data  for  the  case  of  nonholographlc  multiple  glints 
located  on  the  object'.  The  approach  consists  of  three  algorithms 
employed  successively.  The  first  algorithm  reconstructs  the  glints 
only,  both  their  positions  and  their*  complex  values.  It  Involves  the 
triple  Intersection  of  translates  of  the  autocorrelation  function.  The 
second  algorithm  uses  the  reconstructed  glints  along  with  the  aperture- 
plane  Intensity  data  to  arrive  at  a  partially  reconstructed  coherent 
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image  of  the  entire  object.  It  is  modification  of  the  AF-synthes Is 
algorithm  in  x-ray  crystallography.  The  third  algorithm  completes  the 
reconstruction  of  a  high-fidelity  image  using  information  about  the 
support  or  size  of  the  object.  It  is  the  same  iterative  Fourier 
transform  algorithm  that  previously  had  only  been  effective  for  special 
classes  of  coherent  objects;  however,  it  is  effective  for  general 
objects  having  glints  that  are  partially  reconstructed  by  the  first  two 
algorithms. 

In  the  subsections  that  follow,  we  briefly  describe  each  of  the 
three  algorithms  that  make  up  the  imaging  approach  and  show 
reconstruction  results. 

5.2.1  Reconstruction  of  Glints 

In  what  follows  we  assume  that  the.  object  consists  of  both  a  glint 
(or  multiple  glints)  component  g(x),  and  a  diffuse  extended  component, 
d(x) , 

f(x)  ■  g(x)  +  d(x>  (5-1) 

where  x  is  a  2-D  coordinate.  It  Is  assumed  that  the  Fourier  intensity, 

I F(u) 1 2 ,  is  detected,  where  F(u),  G(u)  and  D(u)  are  the  Fourier 

transforms  of  f(x),  g(x)  and  d(x),  respectively.  From  the  Fourier 
Intensity  we  can  compute  the  object's  autocorrelation 

rf(x)  -  [g (x)  ♦  d(x)]  «  [g(x)  +  d(x)]  ■  r’[IF(u)lZ] 

-  g(x)  9  g(x)  +  d(x)  9  d(x)  +  g(x)  9  d(x)  +  d(x)  9  g(x) 

-  rg(x)  +  rd(x)  +  g(x)  9  d(x)  +  d(x)  9  g(«x)  (5-2) 

where  7  denotes  Fourier  transformation,  9  denotes  cross-correlation, 
and  r^  denotes  the  autocorrelation  of  f. 
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If  the  glint  energy  Is 
diffuse  component,  l.e.  If 


K  ■  - w  (5-3) 

X  ld(x)r 

x 

Is  on  the  order  of  one  or  greater,  than  the  peaks  of  rfl(x),  the 
autocorrelation  of  the  glints,  Mill  exceed  the  other  terms  In  Eq. 
(5-2),  enabling  the  glint  Information  to  be  Isolated  from  the  other 
terms  by  a  thresholding  operation  (set  all  values  below  some  threshold 
to  zero),  Once  rg(x)  Is  Isolated,  then  the  glints  can  be  reconstructed 
as  described  below. 

We  model  the  glints  by 

(5-4) 

where  M  Is  the  number  of  glints  and  b(x)  Is  the  Impulse  response  of  the 
Imaging  system  [or  Is  a  delta  function  If  g(x)  represents  the  object]. 
The  autocorrelation  of  the  glints  Is  given  by 

rg(x)  •  g(x)  •  g(x) 

Wn  b(x  '  xm  *  xn>  <5‘5> 

If  the  glints  are  spaced  nonredundantly,  that  Is  If  no  two  vector 
separations  between  distinct  pairs  of  glints  Is  the  same,  then 

rg<xJ  *  «k>  =  Vk  b(0)  *  Vk-  <5'« 


g(x) 


£ 

m«I 


'm 


b(«  -  x„) 


large  compared  with  thfe  energy  of  the 

Z  ig(x)i2 
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where  for  mathematical  simplicity  we  have  normalized  so  that  b(O)  -  1. 
Eq.  (5-6),  which  embodies  the  fact  that  each  glint  In  the 
autocorrelation  arises  from  the  product  of  two  glint  values  In  the 
object,  allows  the  autocorrelation  of  the  glints  to  be  unraveled.  This 
unraveling  can  be  done  In  a  number  of  ways,  as  described  in  Reference 
5.6,  Sections  5  to  7.  The  method  we  used  Is  a  modification  of  one  of 
the  methods  to  which  Ref.  5.6  alludes.  It  consists  of  the  following 
steps. 

1.  To  find  the  glints  In  rf(x)  that  constitute  rg(x),  threshold 
I rf (x) 1 2 3  to  define  an  autocorrelation  support  function  for  r  (x): 


rs(x)  - 


1,  where  I rf(x)lz  i  threshold 


f’ 

lo, 


otherwise 


(5-7) 


Then  we  assume  that  rg(x)  ~  rs(x) rf (x)  except  at  x-0  where  rd(x) 
corrupts  it. 

2.  Reconstruct  the  support  of  the  object  glints  using  the 
autocorrelation  support  trl -intersection  [5.6] 


rsl(x)  ■  r5(x)  rs(x  -  x^,)  r,(x  .  ,.) 


(5-8) 


where  xmaxl  Is  the  location  of  the  maximum  of  I r  (x) I  outside  x-0, 

and  xmaxk  1s  the  1ocation  of  the  maximum  of  lrg(x) I rs(x)r$(x-xmaxl) 
outside  x-0  and 

3.  Multiply  r$j(x)  by  rg(x),  which  can  be  shown  to  yield  either 


rsl(x>  rg(x)  ■  T.  UJ2  b(x)  .  g;  b(x  -  x„  .  xk)  (5- 


9.) 
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■  £  I./  b(x)  ok  2^  «*  b(*  -  *k  +  X„)  (5-9b) 


4.  Using  a  second  support  function,  determine  Ig^l,  which  is  the  value 
at  x-0,  by  the  method  described  in  Section  5.1.5.  This  procedure 
is  required  because  rg(x)  Is  corrupted  at  x*0  by  rd(x). 


5.  Divide  Eq.  (5-9)  by  gk  ■  lgkl  to  get 

glint  (mg.  -  gjj  gy  b(x  -  x,  +  «k) 

(5-10a) 

or 

-  ^  flj  b(x  -  xk  +  xn) 

(5-10b) 

where  (5-10a)  represents  an  image  of  the  glints  (except  gk,  which 
was  already  reconstructed)  and  (5-10b)  represents  a  twin  (complex 
conjugated  and  rotatad  180°)  Image  of  the  glints.  Either  Is  an 
admissible  solution  when  only  the  autocorrelation  or  the  Fourier 
intensity  Is  given. 

This  process  of  reconstructing  the  glints  is  Illustrated  In  Figure 
5-1.  The  object,  In  Fig.  5-1 (a)  Is  a  complex-valued,  speckled  image  of 
a  model  of  a  P72-2  satellite  (the  diffuse  part)  with  three  delta- 
function  glints  artificially  added  on  the  middle  part  of  the  body  of 
the  satellite.  Previously  this  was  thought  to  be  the  most  difficult 
case  for  imaging  using  glints.  (For  a  holographic  reconstruction  there 
must  be  only  one  glint,  and  it  must  be  separated  from  the  body  of  the 
satellite.)  In  this  and  all  the  figures,  the  modulus  of  the  complex- 
valued  Images  are  shown.  Figure  5-l(b)  shows  the  autocorrelation 
function  computed  from  noise-free  Fourier  intensity  data.  Figure 
5-1 (c)  shows  the  results  of  thresholding  the  autocorrelation  function 
Tf(x)  to  isolate  r^(x),  the  autocorrelation  of  the  glints.  M  ■  3 
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FIGURE  5-1.  RECONSTRUCTION  OF  THE  GLINT  COMPONENT  OF  THE  OBJECT.  (A) 
The  object  having  3  glints;  (B)  its  autocorrelation;  (C)  thresholded 
autocorrelation;  (4)  glint  component  of  the  image  reconstructed  from 
(B)  and  (C). 
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glints  in  the  object  produces  MZ-M+1  *  7  glints  In  r*g(x).  The  five- 
step  procedure  described  above  was  then  used  to  successfully 
reconstruct  the  glints  alone,  shown  In  Figure  5-1 (d). 

This  experiment  was  repeated  for  several  different  values  of  the  KJ 
ratio,  as  defined  by  Eq.  (5-3).  The  normalized  root-mean-squared 
(NRMS)  error  of  the  glint  reconstructions  is  plotted  in  Figure  5-2  for 
two  different  cases:  the  three  glints  positioned  on  the  body  of  the 
satellite  as  described  above,  and  the  same  three  glints  positioned  just 
off  the  body  of  the  satellite.  The  positions  of  the  glints  with 
respect  to  the  satellite  body  madr  little  difference.  The  glints  were 
reconstructed  accurately  (NRMS  error  <  0.20)  for  K  *  1.0. 

The  examples  described  above  were  obtained  with  a  simulated  object 
using  simplified  delta-function  glints.  For  the  case  of  realistic 
diffraction-limited  glints  (being  impulse  responses  spread  over  more 
than  one  pixel),  a  modification  of  the  first  step  of  the  5-step 
procedure  Is  required.  It  involves  finding  the  glints  In  the 
autocorrelation  by  a  process  related  to  the  GLEAN  method  from  radio 
astronomy  rather  than  by  simple  thresholding. 

At  this  point  there  are  several  possible  approaches  to  using  the 
reconstructed  glint  Information,  g(x),  to  help  to  reconstruct  the 
entire  object,  f (x) .  In  what  follows  we  describe  only  one  of  the 
approaches,  the  one  that  appeared  to  be  most  effective:  first  obtain  a 
partial  reconstruction  by  AF-synthesIs,  then  complete  the 
reconstruction  by  the  Iterative  transform  algorithm. 

5.2.2  AF  Synthesis  Algorithm 

AF  Synthesis  is  a  method  taken  from  x-ray  crystallographic 
reconstruction,  adapted  by  Baldwin  and  Warner  for  interferometric 
astronomical  Imaging  [5.7]  and  further  adapted  here  for  coherent 
reconstruction  when  glints  are  known.  It  works  poorly  when  little 
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about  the  Image  Is  known,  but  works  well  when  parts  of  the  Image  are 
known  well.  This  Is  the  case  If  the  glints  are  reconstructed  as 
described  In  the  previous  section. 

XL 

One  Iteration  of  AF  synthesis  Is  as  follows.  At  the  n  Iteration 
we  assume  we  know  some  part  gn(x),  of  the  Image,  where 

f(x)  -  gn(x)  +  dn(x)  (5-11) 


and  dn(x)  Is  the  unknown  part  of  the  Image.  For  the  first  Iteration 
g^x)  Is  the  Image  of  the  glints  reconstructed  by  the  method  described 
In  the  previous  section  and  dj(x)  Is  the  (unknown)  diffuse  part  of  the 
object.  The  modulus  squared  of  the  Fourier  transform  of  Eq.  (5-11)  Is 


IF(u)l2  ■  IGn(u)l2  +  gJ(u)  Dn(u)  +  Gn(u)  D*(u)  +  IDn(u)l2  .  (5-12) 

Taking  the  Fourier  transform,  Gn(u),  of  gn(x)  and  multiplying  It  by 
I F (u) I / I Gn (u) I  yields 

& 

Gn  '  Gn  'F'/"V 

’  K'*  *  «X  +  5n°n  +  V'V  . 

=  G„  ♦  0/2  *  (S„  'V2  +  GnDX2IGn|2)  •  <►«) 

using  a  Taylor-serles  expansion  assuming  IGnl  large.  Therefore 


2<  •  Gn  =  sn  +  °n  *  (S  10,2  +  GM/I5n'2  ■ 

. . .  "I 


(8-14) 
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Inverse  Fourier  transforming  yields 

2g^(x)  -  gn(x)  -  f(x)  +  other  terms.  (5-15) 

That  Is,  by  subtracting  the  Input  Image,  gn(x),  from  twice  the  output 
image,  g^(x),  we  reconstruct  the  entire  image,  f(x).  However,  the 
"other  terms"  severely  corrupt  much  of  the  desired  image. 
Nevertheless,  the  brightest  points  in  Eq.  (5-15)  outside  the  previously 
known  image  points,  gn(x),  are  most  likely  to  belong  to  the  object, 
f(x),  rather  than  to  the  other  terms.  Thus  we  take  the  brightest  new 
points  resulting  from  the  computation  of  Eq.  (5-15)  and  add  those  to 
gn(x)  to  form  gn+1(x)  which  represents  a  larger  known  portion  of  the 
image.  This  is  done  repeatedly  until  a  reasonably  complete  partially- 
reconstructed  image  appears. 

Figure  5-3  shows  the  results  of  the  application  of  AF  synthesis. 
Figure  5-3 (a)  shows  the  same  coherent  object  with  three  glints  on  the 
body,  and  Figure  5-3 (b)  shows  the  same  reconstructed  glints  as  for 

i 

Figure  5-1.  (Any  difference  in  appearance  from  Figure  5-1  is  due  to  a 
difference  in  exposures  when  photographing  the  images.)  Using  the 
reconstructed  glints  as  g^x),  Eq.  (5-15)  was  computed,  the  result  of 
which  is  shown  in  Figure  5-3(c).  Much  of  the  object  is  apparent  In 
Figure  5-3 (c) ,  but  there  is  a  high  level  of  artifacts.  Figure  5-3 (d) 
shows  the  improved  results  after  12  iterations  of  AF  synthesis.  The 
locations  of  the  points  assumed  to  be  known  at  this  stage  are  shown  in 
Figure  5-3(e).  Further  iterations  of  AF  synthesis  resulted  in  little 
further  improvement. 

5.2.3  Iterative  Fourier  Transform  Algorithm 

The  iterative  Fourier  transform  algorithm  involves  repeatedly 
Fourier  transforming  an  estimate  back  and  forth  between  the  Fourier 
domain,  where  the  measured  Fourier  intensity  is  reinforced,  and  the 
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FIGURE  5-3.  IMAGE  RECONSTRUCTION  BY  THE  THREE-ALGORITHM  METHOD.  (A) 
The  object;  (B)  the  reconstructed  glints  (from  Figure  5-1);  (C)  output 
from  one  iteration  of  AF  synthesis;  (D)  output  from  12  iterations  of  AF 
synthesis;  (E)  locations  of  points  assumed  to  be  known  in  (D);  (F) 
image  reconstructed  by  the  iterative  transform  algorithm  staring  with 
(D). 
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image  domain,  where  a  support  constraint  is  enforced.  The  support 
constraint,  i.e.  knowledge  that  the  image  is  zero  outside  some  well- 
defined  area,  can  be  derived  from  the  autocorrelation  of  the  object 
[5.6].  The  details  of  the  Iterative  Fourier  transform  algorithm  are 
described  elsewhere  [5.2, 5. 8, 5. 9] .  The  partially  reconstructed  Image 
shown  in  Figure  5-3 (d)  was  used  as  an  Input  to  the  Iterative  Fourier 
transform  algorithm  for  which  we  employed  a  supports  constraint  in  the 
form  of  a  crude  rectangle  Inside  of  which  the  object  loosely  fit.  The 
image  reconstructed  by  the  Iterative  Fourier  transform  algorithm  Is 
shown  In  Figure  5-3 (f) .  It  Is  an  excellent  reconstruction  of  the 
object  shown  in  Figure  5-3 (a),  despite  having  assumed  total  loss  of 
phase  Information  and  despite  the  fact  that  the  glints  were  In  the  most 
unfavorable  locations.  (Prior  to  this  work  It  was  thought  that 
multiple  glints  Imbedded  within  the  object  would  be  the  most  difficult 
case.) 

This  experiment  was  repeated  for  several  K-ratlos  (glint  energies), 
and  the  results  are  shown  In  Figure  5-4  and  plotted  In  Figure  5-5. 
Very  recognizable  images  were  reconstructed  for  K  i  2  (Image  NRMS  error 
i  0.35).  For  K  •  5,  the  simulation  experiment  was  repeated  for  various 
levels  of  photon  noise  In  the  Intensity  data.  The  results  are  shown  In 
Figure  5-6  and  plotted  In  Figure  5-7.  Recognizable  Images  were 
reconstructed  for  107  or  greater  photons  per  intensity  array  of  128  x 
128  samples  (or  about  60  x  40  «  3200  speckles),  equivalent  to  about  600 
photons  per  sample  or  4000  photons  per  speckle. 

5.2.4  Summary  of  the  Three-Algorithm  Method 

Prior  to  this  work,  coherent  image  reconstruction  from  a  single 
snapshot  of  far-field  laser  speckle  intensity  data  was  possible  if  the 
object  included  a  single  well -separated  (holographic)  glint  or  a  very 
bright  glint  beyond  the  edge  of  the  diffuse  part  of  the  object.  We 
have  developed  an  approach  for  reconstructing  objects  having  much  less 
favorable  glints,  including  multiple  glints  that  may  be  located  within 
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FIGURE  5-4.  IMAGES  RECONSTRUCTED  BY  THE  THREE-ALGORITHM  METHOD  FOR 
VARIOUS  GLINT  ENERGIES  (K-RATIOS).  (A) -(D)  Three  glints  off  the 
object,  (E)-(H)  three  glints  on  the  object. 
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FIGURE  5-5.  NORMALIZED  RMS  ERROR  IN  RECONSTRUCTING  THE  IMAGE  AS  A 
FUNCTION  OF  K  RATIO. 
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FIGURE  5-6.  IMAGES  RECONSTRUCTED  BY  THE  THREE-ALGORITHM  METHOD  FOR 
VARIOUS  PHOTON  LEVELS.  K“5  and  three  glints  on  object. 
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FIGURE  5-7.  NORMALIZED  RMS  ERROR  VS.  NUMBER  OF  PHOTONS. 


the  diffuse  part  of  the  object.  The  approach  consists  of  three 
successive, algorithms:  (1)  a  triple  intersection  of  the  autocorrelation 
function  that  yields  an  image  of  the  glints  alone,  (2)  the  AF-synthesis 
algorithm  that  uses  the  image  of  the  glints  together  with  the  Fourier 
intensity  data  to  yield  a  partial  reconstruction  of  the  entire  image, 
and  (3)  the  iterative  Fourier  transform  algorithm  that  uses  the 
partially  reconstructed  image  together  with  the  Fourier  Intensity  data 
and  a  support  constraint  to  complete  the  reconstruction.  For  the 
example  investigated  having  three  delta-function  glints,  good 
reconstructions  were  obtained  for  K  -  (glint  energy) /(diffuse/energy)  i 
2  and  4000  photons  per  detected  speckle  (i.e.,  a  relatively  high  light 
level). 

Further  research  is  required  to  optimize  the  approach  and  to 
quantify  performance  for  diffraction-limited  (as  opposed  to  delta- 
function)  glints,  to  extend  the  method  to  work  for  extended  glints,  and 
to  demonstrate  the  method  on  laboratory  experimental  and  field  data. 

5.2.5  Determining  Ig^l  for  Glint  Reconstruction 

To  determine  the  value  of  the  first  glint,  Ig^l,  perform  the 
following  steps. 

1.  Sum  over  the  squared  modulus  of  rsi(x)**g(x)  for  m^k: 

ck  ■  's/  'a/  ■  <z:  '<>/  ■  '■/  •  (5-16> 

2.  Reconstruct  a  second  image  support: 

rs2(x)  -  rs(x>  rs(x  -  x^,)  r,(xraxj)  (5-17) 

where  x[naxj  is  the  location  of  the  largest  peak  of  r$(x)rs(x  - 
xmaxi)rg(x)  that  is  outside  r$j(x). 
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3.  As  in  Step  1,  sum  over  the  squared  modulus  of  r^r^x)  for  ni*j: 

Cj  ■  I8ji2  2~  ig„i2  ■  i9ji2nZ  ig„i2  -  Is/  •  (s-is) 


4.  Let 


cjk  •  r<W  ■  r<*J  '  xk>  ’  Ojflk 


(5-19) 


5.  Solve  Eqs,  (5-16),  (5-18)  and  (5-19)  for  lgklz: 


lgklz  -  iCjkl 


CL  -  1C 


Ik1 


LCj  -  'Cjk1' 


1/2 


6.  Assume  that  gk  has  zero  phase  (we  can  arbitrarily  set  any  one  phase 
value  to  whatever  we  want): 


gk  -  ilgkl2  .  (5-20) 

Note  that  a  simple  method  Is  possible  If  the  object  had  no  diffuse 
part;  however,  the  diffuse  part  of  the  object  will  usually  make  a 
strong  contribution  to  r^(0),  making  this  5-step  method  necessary. 

5.3  RECONSTRUCTION  WITH  ONLY  THE  ITERATIVE  TRANSFORM  ALGORITHM 

Although  the  three-algorithm  method  discussed  In  th?  previous 
section  Is  generally  the  most  robust  way  to  reconstruct  an  object 
having  glints  from  a  single  snapshot  of  Fourier  intensity  data,  It  Is 
also  possible  to  reconstruct  using  only  the  Iterative  transform 
algorithm  (the  third  algorithm).  In  this  case  It  Is  assumed  that  by 
Inspection  of  the  autocorrelation  function,  computed  from  the  measured 
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Fourier  Intensity,  we  can  determine  that  there  are  one  or  more  glints 
in  the  object.  Furthermore,  If  there  are  two  glints,  then  the 
separation  of  the  glints  can  also  be  determined  from  the 
autocorrelation. 

To  start  the  Iterative  transform  algorithm,  two  different  starting 
estimates  were  tried.  First,  since  the  autocorrelation  function 
contains  the  desired  image  as  one  of  its  terms  (as  described  In  Sect. 
5.3),  we  chose  as  an  initial  estimate  a  windowed  version  of  the 
autocorrelation  function.  Second,  we  chose  the  initial  estimate  to  be 
a  bright  glint  (or  glints)  surrounded  by  random  noise.  In  practice  the 
latter  initial  estimate  worked  better.  It  was  assumed  that  the 

approximate  location  of  the  glint  (or  glints)  was  known. 

.1 

Figure  5-8  shows  reconstruction  results  for  the  case  of  a  single 
bright  glint  added  near  the  top  center  of  the  speckled  satellite  image. 
The  support  constraint,  assumed  known  a  priori,  was  a  rectangle' just 
enclosing  the  object  with  the  glint.  A  plot  of  the  nonhalized  RMS 
error  of  the  reconstructed  Images  as  a  function  of  the  K  ratio  is  shown 
in  Figure  5-9.  The  reconstructions  are  very  good  for  K  greater  than 
0.25.  These  reconstructions  are  better  than  those  of  the  previous 
section  because  (1)  since  there  is  only  one  glint,  that  single  glint 
contains  all  the  energy  implied  by  the  K  ratio,  making  it  relatively 
brighter  than  any  of  the  three  glints  for  the  previous  case,  for  any 
given  value  of  K,  and  (2)  reconstruction  from  three'  glints  is 
inherently  much  more  difficult  than  from  a  single  glint.  For  the  cases 
of  K  ■  1.0  and  0.5,  Poisson  noise  was  added  to  the  Fourier  intensity 
data.  The  errors  in  the  corresponding  reconstructed  Images  is  shown  in 
Figures  5-10(a)  and  (b)  respectively.  ■ 

Figure  5-11  shows  reconstruction  results  for  the  case  of  a  pair  of 
glints,  one  of  which  is  in  the  same  position  as  the  case  above,  and  the 
second  is  nearer  to  the  center  of  the  object.  Figure  5-12  shows  a  plot 
of  the  normalized  RMS  error  as  a  function  of  glint  strength.  Very  good 
reconstructions  are  obtained  for  K  >  0.5. 
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FIGURE  5-8.  IMAGES  WITH  A  SINGLE  GLINT  RECONSTRUCTED  BY  THE  ITERATIVE 
TRANSFORM  ALGORITHM,  FOR  VARIOUS  K  RATIOS. 
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FIGURE  5-10.  NORMALIZED  RMS 
SINGLE  GLINT  RECONSTRUCTED  BY 
For  K*1.0;  (b)  for  K-0.5. 


ERROR  VS.  DATA  ERROR  FOR  IMAGES  WITH  A 
THE  ITERATIVE  TRANSFORM  ALGORITHM,  (a) 


FIGURE  5-11.  IMAGES  WITH  TWO  GLINTS  RECONSTRUCTED  BY  THE  ITERATIVE 
TRANSFORM  ALGORITHM,  FOR  VARIOUS  K  RATIOS. 
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FIGURE  5-12.  NORMALIZED  RMS  ERROR  VS.  K  RATIO  FOR  IMAGES  WITH  TWO 
GLINTS  RECONSTRUCTED  BY  THE  ITERATIVE  TRANSFORM  ALGORITHM. 


Diffraction  effects  due  to  the  aperture  are  not  included  in  this 
simulation.  The  inclusion  of  diffraction  effects  would  yield  poorer 
performance;  further  research  is  needed  to  develop  improved  algorithms 
for  this  case. 

5.4  RECURSIVE  AUTOCORRELATION  ALGORITHM 

Initially  a  recursive  autocorrelation-based  reconstruction 
algorithm  was  investigated.  It  was  later  overshadowed  by  the 
algorithms  described  above,  but  we  include  it  here  for  the  sake  of 
completeness. 

We  assume  that  the  object  has  a  single  glint  of  magnitude  a  and  a 
diffuse  part.  d(x),  so  that  it  can  be  modeled  as 

f(x)  ■  a  d(x)  +  d(x)  (5-21) 

and  its  autocorrelation  function  as 

rf(x)  ■  I al 2  d(x)  +  rd(x)  +  a*  d(x)  +  a  d*(-x)  (5-22) 

which  contains  within  it  a  representation  of  the  diffuse  part,  d(x), 
multiplied  by  a*.  Unless  the  glint  satisfies  the  holography  condition, 
however,  the  other  terms  will  overlap  the  desired  term,  a*  d(x),  making 
It  not  Immediately  available.  If  the  glint  does  not  satisfy  the 
holography  condition  but  Is  to  one  side  of  d(x),  then  only  the  term 
rd(x)  overlaps  the  desired  term.  Under  that  circumstance,  the 
following  recursive  reconstruction  algorithm  is  possible.  After 
estimating  the  glint  strength,  first  estimate  a*  d(x)  by  windowing  one 
side  of  r^(x).  Then  form  a  second  estimate  by  subtracting  from  r^(x) 
an  estimate  of  rd(x),  computed  from  the  estimate  of  d(x),  and  windowing 
the  result.  This  process  Is  continued  until  no  further  changes  are 
made.  In  practice  this  method  requires  a  strong  glint  and  was  not  as 
successful  as  the  methods  described  ip  the  previous  sections. 
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5.5  EFFECT  OF  GLINT  STRENGTH  ON  DATA  QUANTIZATION  ERROR 

It  Is  assumed  that  the  laser-illuminated  object  has  a  glint 
component  and  a  diffuse  component.  It  may  be  that  the  glint,  say  from 
a  flat  panel,  appears  only  for  certain  angular  orientations  of  the 
object.  A  problem  is  that  the  energy  from  the  glint  can  far  exceed  the 
energy  from  the  entire  diffuse  component  of  the  object.  Then  aperture- 
plane  detectors  having  a  finite  dynamic  range  (or  a  finite  number  of 
quantization  levels)  may  have  the  Information  about  the  diffuse 
component  overwhelmed  by  the  energy  from  the  glint.  In  this  section 
this  problem  Is  analyzed. 

For  the  case  of  a  single  glint,  the  model  for  the  object  Is  again 

f(x)  •  a  5(x)  +  d(x)  (5-23) 

which  has  Fourier  transform 


F(u)  -  a  +  D(u)  (5-24) 

2 

where  D(u)  Is  a  zero-mean  Gaussian  random  variable  with  variance  a. 
The  detected  quantity  Is  the  Intensity 

I F(u) I2  -  la  +  D(u) 1 2 

■  a2  +  2  a  Re[D(u)]  +  I 0(u) I2  ,  (5-25) 

where  we  have  used  the  fact  that  the  phase  of  the  glint,  a,  can, 
without  loss  of  generality,  be  set  to  zero.  Letting  w  *  I F (u) 1 2  and  g 
■  a2,  the  intensity  follows  a  modified  Rician  or  non-central  Wishart 
distribution: 

pw(w)  *  (2*)"1  exp[-w/(2e2)]  exp [-g/ (2a2)]  I0(W/ff2)  •  (5-26) 
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To  simplify  this  further,  we  can  define  intensities  normalised  to  the 

energy  In  the  diffuse  part  of  the  objects  replacing  w/(2 r)  by  w  and 

2  •  .  '  *  '  #  ’ 

g/(2 a  )  by  g',  we  have 

.  i  1  , 

Pw(w)  •  e"w  e’fl  I0(2JwT)  . ,  (5-27) 

The  mean  and  variance  of  w  are  1  +  g  and  1  +  2  g,  respectively. 

t  i  ‘  ■  ■  1 

,  .  .  i  -  -»  •  i 

Jr  1  * 

,  l  ' 

The  assumed  quantization  operation  Is  Illustrated  In  Figure  5-13, 
which  shows  a  linear  quantizer  with  an  offset  7  and  N  quantization 
Intervals,  each  of  width  Aw.  For  a  given  value  of  the  relative  glint 
strength,  g,  the  mean-squared  error  due  to  the  quantization  rule  on  the 
probability  distribution  pw(w)  can  be  numerically  computed. 

Optimum  quantizers  were  derived  for  two  casess  no  glint  and  a  very 
strong  glint.  This  was  accomplished  by  arriving  at  a  slgnal-to-nolse  , 
expression  as  a  function  of  N,  Aw,*  and  7,  and  then  solving  nonlinear 
equations  for  the  optimum  Aw  and  7  using  Newton's  method.  The  results 
are  summarized  In  Table  5-1. 


Optimum  Quantizer  Intervals 

For  zero  allnt 

.  '  1, 

For  larae  allnt 

N 

r 

SRN(db) 

_ Aw 

SRN(db) 

256 

.04368 

37.27 

.03076 

40.57 

512 

.02432 

42.43 

.01650 

46.04 

1024 

.01341  . 

47.65 

.00879 

51.55 

2048 

.00734 

52*94 

. .00465 

67.11 

4096 

.00399 

58.28 

.00245 

62.71 

The  slgnal-to-nolse  ratios  (SNR's)  of  four  quantizers  were  compared 

I  I  >  I  V 

for  the  case  of  a, design  optimized  for  g  *  10 ,0d0  (an  enormous  maximum 
glint  strength),  for  various  values  of  N  from  256  to  4096,  and  for 
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5-13.  UNIFORM  QUANTIZER  FOR  MODIFIED  RICIAN  LASER 


values  of  an  actual  glint  strength  of  0,  10,  100,  1,000,  and  10,000. 
The  results  are  shown  In  Table  5.2.  From  this  table  it  can  be  seen 
that  with  a  fixed  step  size  (optimized  for  g  ■  10,000),  the  SNR  for  the 
smaller  values  of  g  Is  very  poor.  The  addition  of  a  dynamic,  optimum 
step  size  yields  dramatic  Improvement  In  SNR  for  the  smaller-gllnt 
cases.  The  same  effect  as  a  dynamic  step  size  can  be  achieved  by  an 
automatic  gain  control  (say,  attenuating  the  llfcht  arriving  at  the 
detector  when  the  glint  becomes  very  bright).  The  addition  a  dynamic, 
optimum  offset  yields  additional  gains  In  SNR,  especially  for  the 
brlghter-gilnt  cases. 

The  reason  that  automatic  gain  control  Is  Important  Is  seen  from 
the  fact  that  In  the  second  term,  the  Information-carrying  quantity,  of 
Eq.  (5-15)  Is  multiplied  by  the  glint  strength. 

The  major  conclusion  from  this  study  Is  that  In  order  to  allow  for 
very  large  glints,  we  should  have  a  dynamic  or  adaptive  quantizer.  The 
most  Important  feature  of  the  dynamic  quantizer  would  be  a  variable 
step  size  Aw.  Thl?  could  be  achieved  by  an  automatic  gain  control. 
Helpful,  but  less  Important,  Is  to  allow  for  a  dynamic  offset  7.  The 
results  were  derived  for  very  large  glints;  the  optimization  of  the 
quantizer  for  Intermediate-strength  glints,  which  are  likely  to  occur 
In  practice,  would  require  additional  research. 
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Table  5-2.  Comparison  of  Uniform  Quantizers.  The  tables  give  the  SNR  as  a  function  of  relative 
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6.0  IMAGING  WITH  PARTIAL  PHASE  INFORMATION 


In  this  section  and  in  Section  7,  we  describe  methods  developed  for 
using  partial  phase  information  in  the  phase  retrieval/image 
reconstruction  process.  For  the  most  difficult  objects  to  reconstruct 
(complex-valued,  having  no  glints  or  separated  parts),  some  additional 
information  is  essential  to  obtain  a  reliable  reconstruction,  given  our 
present  algorithms.  One  kind  of  such  additional  information  is  partial 
phase  information.  Partial  phase  information  can  pome,  for  example, 
from  an  Imaging  system  that  Inherently  measures  or  determines  the 
phase.  If  the  partial  phase  Information  is  nearly  complete,  so  that 
using  that  phase  yields  a  useful  image,  then  the  phase-retrieval 
processing  can  be  thought  of  as  a  way  t:  clean  up  the  image  to  Improve 
its  quality.  This  is  equivalent  to  reducing  the  errors  or  filling  In 
the  gaps  In  the  given  partial  phase  information.  If  the  partial  phase 
information  is  very  incomplete  or  noisy,  then  no  useful  image  would 
result  from  it,  and  the  phase-retrieval  processing  would  be  forming  the 
Image.  In  the  first  place,  with  the  partial  phase  information  helping  it 
to  succeed. 

Two  major  cases  of  partial  phase  information  were  considered:  (1) 
phase  known  well  over  a  small  aperture,  and  (2)  noisy  phase  over  the 
entire  aperture.  In  the  first  case,  it  is  a  matter  of  filling  In  the 
missing  phase,  but  most  of  the  phase  Is  missing.  In  the  second  case, 
it  Is  a  matter  of  correcting  the  errors  in  the  given  phase. 

Two  scenarios  that  would  correspond  to  the  first  case  are  as 
follows.  Suppose  that  the  object  Is  coherently  Illuminated  with  a 
laser,  and  Intensity  measurements  are  made  In  an  aperture  plane  of  the 
optical  system.  In  addition,  optical  field  measurements  are  made  over 
a  smaller  aperture  imbedded  In  (or  contiguous  with)  the  Intensity 
measurements.  The  optical  field  measurements  could  be  performed,  for 
example,  by  heterodyne  detection  or  by  two  Intensity  measurements  In 
different  planes  and  the  fields  reconstructed  by  the  Gerchberg-Saxton 
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algorithm.  These  optical  fields  would  be  aberrated  by  atmospheric 
turbulence.  In  addition,  the  small  aperture  would  also  have  a 
wavefront  sensor  that  measures  the  atmospheric  wavefront  error. 
Aberration-free  optical  field  data  could  be  obtained  over  the  small 
aperture  by  subtracting  the  phase  due  to  the  atmosphere  from  the  phase 
of  the  measured  optical  field. 

The  secohd  scenario  representing  the  first  case  would  be  for  a 
system  in  space  having  a  small,  diffraction-limited  telescope  making 
optical  field  measurements,  Imbedded  In  a  larger  aperture  over  which 
the  aperture-plane  intensity  measurements  are  made.  In  this  case  a 
wavefront  sensor  would  not  be  needed  since  there  would  be  no 
atmospheric  turbulence  to  aberrate  the  optical  data  for  the  small, 
diffraction-limited  telescope. 

Many  different  imaging  systems  could  provide  data  for  the  second 
case,  that  of  noisy  phase  known  ovter  the  entire  aperture.  They  include 
active  imaging  'modalities  such  as  the  Itek/Western  system,  triple 
correlation  of  aperture-plane  Intensity,  and  FOCI  and  passive  imaging 
modalities  such  as  astronomical  speckle  interferometry  using  triple 
correlation  and  aperture-plane  Interferometry  using  phase  closure. 

The  basic  approach  to  phase  retrieval  and  Image  reconstruction 
taken  for  these  scenarios  was  to  use  the  iterative  Fourier  transform 
algorithm  to  take  advantage  of  all  the  available  data  and  constraints 
to  form  the  solution.  This  approach  allows  for  the  combination  of  a 
variety  of  disparate  types  of  Information,  such  as  Fourier  modulus 
(square  root  of  intensity),  Fourier  phase,  object-domain  support 
(finite  extent)  constraint,  and  nonnegattvlty  (applicable  for 
Incoherent  Images).  For  the  case  of  phase  know  well  over  a  small 
aperture,  a  variation  of  the  Iterative  transform  algorithm,  called  the 
expanding-weighted-mask  algorithm,  was  developed.  This  Initial  attempt 
with  the  expand!  ng<-weighted-mask  algorithm  jjave  results  that  were  only 
partially  successful:  they  were  promising  but  very  preliminary  and 
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Incomplete.  (In  a  separate,  later  program  funded  by  the  Naval  Research 
Laboratory,  this  approach  was  expanded  upon  and  optimized,  and  yielded 
very  good  reconstructed  Images.)  Those  preliminary  results  are 
described  in  Section  6.1.  For  the  case  of  noisy  phase  over  the  entire 
aperture,  a  variation  of  the  Iterative  transform  algorithm,  called  the 
phase  variance  algorithm,  was  developed,  as  described  in  Section  6.3. 
In  addition,  an  alternative,  entirely  new  approach  was  developed  for 
that  case:  2-D  shear  averaging,  which  Is  described  In  Section  7.  One 
other  case  that  was  briefly  explored  was  that  of  knowing  one  bit  of 
phase.  In  that  case,  reconstruction  was  easily  achieved  with  a 
windowing  of  the  Initial  Image  computed  from  the  given  phase  followed 
by  cleaning  up  with  the  Iterative  Fourier  transform  algorithm.  This 
last  case,  mostly  of  academic  Interest  since  It  does  not  naturally 
occur  In  currently  known  Imaging  sensors,  Is  described  In  Section  6.2. 

6.1  THE  EXPANDING  WEIGHTED  MODULUS  ALGORITHM 

If  the  Fourier  Intensity  Is  measured  over  a  large  aperture  and  the 
phase  Is  measured  over  a  small  aperture  imbedded  In  the  large  aperture, 
then  It  Is  possible  to  use  that  known  phase  to  help  to  retrieve  the 
phase  over  the  large  aperture.  This  can  be  accomplished  by  enforcing 
the  known  phase  together  with  all  the  other  available  Information 
(Fourier  modulus,  object  support  constraint)  using  the  iterative 
Fourier  transform  algorithm. 

The  support  constraint  can  be  gotten  from  the  available  data  In  one 
of  two  ways.  First,  one  can  use  a  triple  Intersection  of  the 
autocorrelation  support  computed  from  the  Fourier  Intensity  to  put  an 
upper  bound  on  the  support  of  the  object  [6.1].  Second,  from  the 
small -aperture  phase  combined  with  the  measured  intensity  over  the 
small  aperture,  one  can  gets  a  diffraction-limited,  but  low  resolution 
(owing  to  the  small  size  of  the  aperture)  Image.  A  support  constraint 
can  be  formed  by  an  appropriate  thresholding  of  this  low-resolution 
Image. 
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The  Initial  estimate  for  the  Iterative  transform  algorithm  can  be 
gotten  by  simply  using  the  complex-valued  low-resolution  image  or  by 
filling  the  support  constraint  with  complex-valued  random  numbers. 

When  the  Iterative  transform  algorithm  was  run  with  either  of  the 
two  initial  estimates  and  using  either  of  the  two  support  constraints, 
and  enforcing  the  small -aperture  phase,  it  stagnated  without  converging 
to  a  solution.  Essentially  random  phases  were  produced  outside  the 
small  aperture.  Since  the  ratio  of  the  area  of  the  large  aperture  to 
that  of  the  small  aperture  was  chosen  to  be  a  large  number,  the  random 
phases  outside  the  small  aperture  overwhelmed  the  Influence  of  the 
correct  phases  within  the  small  aperture. 

In  order  to  combat  this  problem,  we  began  development  of  the 
expanding  weighted  modulus  algorithm.  It  consists  of  the  following 
steps.  First  the  Fourier  modulus  (the  square  root  of  the  measured 
intensity)  Is  multiplied  (weighted)  by  an  apodlzing  function  that  goes 
to  zero  over  an  area  only  somewhat  larger  than  the  area  of  the  small 
aperture.  Then  several  iterations  are  performed.  The  idea  is  that 
with  the  weighting  function  in  place,  the  known  phase  will  not  be 
overwhelmed  by  the  unknown  phase,  which  now  exists  over  a  much  smaller 
area  than  before,  and  furthermore  has  an  associated  magnitude  that  Is 
weighted  down  in  the  area  of  the  unknown  phase,  further  reducing  Its 
Influence.  Thus  the  known  phase  has  a  chance  to  be  useful  as  a 
constraint  that  helps  to  retrieve  the  unknown  phase  over  the  larger 
area.  Next  the  weighting  function  Is  replaced  by  a  weighting  function 
that  Is  nonzero  over  a  wider  area.  Then  more  iterations  are  performed, 
retrieving  the  unknown  phase  over  this  wider  area.  This  process  of 
widening  the  weighting  function  and  performing  more  iterations  Is 
repeated  until  the  weighting  function  Is  nonzero  over  the  entire  large 
aperture,  at  which  point  the  entire  phase  is  retrieved  and  a  fine- 
resolution  image  is  reconstructed. 
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Several  trials  of  the  expanding  weighted  modulus  algorithm  were 
made,  and,  while  the  results  showed  promise,  the  images  that  were 
produced  were  far  from  diffraction  limited.  For  these  experiments,  the 
weighting  functions  were  chosen  to  be  either  rectangle  functions  or 
triangle  functions,  and  the  number  of  intermediate  weighting  functions 
used  was  small.  (These  initial  results  were  greatly  Improved  upon  In  a 
separate  effort,  In  which  it  was  found  that  by  using  weighting 
functions  with  continuous  derivatives  and  using  a  very  much  larger 
number  of  intermediate  weighting  functions,  the  reconstructions  could 
be  reliable  and  of  high  quality  [6.2]). 

6.2  RECONSTRUCTION  WITH  ONE  BIT  OF  PHASE 

One  bit  of  phase  information,  which  is  equivalent  to  knowing  the 
sign  of  the  real  part  of  the  Fourier  transform,  is  well  known  to 
contain  considerable  information.  First  of  all,  if  the  image  is 
"causal,"  i.e.,  it  Is  located  completely  to  one  side  of  the  optical 
axis,  then  the  inverse  Fourier  transform  of  the  Fourier  modulus 
combined  with  the  one  bit  of  phase  yields  an  image  plane  with  th« 
following  components:  the  desired  image,  the  twin  (complex  conjugated 
and  reflected  about  the  origin)  of  the  desired  Image,  and  noise  and 
artifacts.  The  strength  of  the  noise  and  artifacts  depends  on  the 
degree  of  oversampling  of  the  Fourier  modulus  data.  With  a  high  degree 
of  oversampling  the  noise  and  artifact  level  can  be  very  low,  yielding 
a  good-quality  image  immediately.  This  is  Illustrated  in  Figure  6-1, 
which  shows  the  object  (a)  and  the  image  reconstructed  by  inverse 
transforming  the  Fourier  magnitude  plus  one  bit  of  phase  (c) .  For 
comparison,  Figure  6-1 (b)  shows  the  inverse  transform  of  a  constant 
modulus  with  the  one  bit  of  phase.  The  quality  of  this  latter  image 
shows  that  the  one  bit  of  phase  without  any  modulus  information  is  very 
useful  indeed.  Figure  6-1 (d)  shows  the  result  of  using  a  windowed 
version  of  Figure  6-1 (c)  as  an  initial  estimate,  then  performing 
several  Iterations  of  the  iterative  transform  algorithm  using  a  support 
constraint  In  order  to  refine  the  phase.  The  iterative  transform 
algorithm  successfully  removed  most  of  the  noise  and  artifacts. 
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FIGURE  6-1.  RECONSTRUCTION  WITH  ONE  BIT  OF  PHASE,  (a)  Object,  (b) 
image  from  a  constant  modulus  combined  with  one  bit  of  phase,  (c)  image 
from  the  correct  modulus  combined  with  one  bit  of  phase,  (d)  image 
reconstructed  by  the  iterative  Fourier  transform  algorithm  using  (c)  as 
an  initial  estimate. 


From  Figure  6-1  we  see  that  one  bit  of  phase  Information  Is  very 
powerful  Information,  and  that  what  noise  and  artifacts  It  Introduces 
can  be  easily  cleaned  up  by  the  Iterative  Fourier  transform  algorithm. 
Unfortunately,  none  of  the  sensors  presently  under  development  yield  an 
accurate  measure  of  one  bit  of  phase,  so  the  results  shown  In  Figure 
6-1  are  presently  of  academic  Interest  only. 

6.3  PHASE  VARIANCE  ALGORITHM 

In  what  follows  Is  described  a  modification  of  the  Iterative 
transform  algorithm  which  uses  a  poorly-known  phase  across  the  Fourler- 
domaln  aperture. 

Let  the  object  and  Its  Fourier  transform  be  f (x)  and  F(u)  ■  IF(u)l 
exp[1*(u)] ,  respectively.  Suppose  we  measure 

G0(u)  -  IG0(u)l  exp[10(u)] 

■  F(u)  exp[1^e(u)]  -  IF(u)l  exp{1[*(u)  +  #e(u)]>  (6-1) 

where  ^  (u)  Is  a  phase  error  with  known  (or  known  approximately) 

v 

variance  a.  So  the  measured  (noisy)  phase  Is 

fi(u)  »  [f (u)  +  *e(u)]mod  2f  •  (6-2) 

The  Image  gotten  by  Inverse  Fourier  transforming  GQ(u)  would  be  g0(x), 
a  blurred  Image. 

We  seek  ways  to  Improve  the  phase  estimate  over  that  given  by  the 
measurement  0(u).  This  may  be  accomplished  by  the  Iterative  Fourier 
transform  phase  retrieval  algorithm  which  uses  additional  Information 
In  the  object  domain,  such  as  nonnegativity  and/or  support  constraints 
to  Infer  the  true  phase  of  F(u).  Two  approaches  are  described  next. 
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The  first  approach  Is  tb  perform  the  usual  phase  retrieval 
algorithm,  typically  cycles  of  hybrid  Input-output  (HID)  and  error- 
reduction  (ER),  and  simply  use  0(u)  as  the  Initial  estimate  for  the 
Fourier-domaln  phase.  The  Fourler-domafn  constraint  would  be  the 

measured  modulus  1 GQ(u) I  »  I F (u) I . 

1“ 

The  second  approach  Is  to  constrain  the  phase  during  the  Iterations 
to  lie  near  6(u).  Constraining  the  phase  to  equal  6( u)  does  no  good 
since  one  would  simply  get  the  blurred  Image  with  no  change.  Instead 
it ■ is  more> useful  to  allow  the  phase  to  wander  from  tf(u),  but  not  let 

.  t  ■ , ' 

It  wander  too  far.  This  can  be  accomplished  using  the  phase  variance 
algorithm,  which  Is  described  as  follows.  In  the  Fourier  domain,  as 
well  as  constraining  the  modulus  to  equal  IF(u)l,  constrain  the  phase 
to  not  differ  from  0(u)  by  more  than  co,  where  c,  the  variance  factor, 
Is  a  real  constant  on  the  order  of  unity.  In  order  to  account  for  2ir 
ambiguities,  this  should  be  performed  as  follows,  where  ^  Is  the  phase 
of  the  Fourier  transform  of  the  input  object  to  the  Iterative  loop  and 
Is  the  altered  phase: 


0  -  co 

•  •  S>mod  2,  < 

r - 1 

* 

'  s  «  ’  9>nmd  2r  s  » 

(6-3) 

0  +  co 

•  V  ‘  2r  >  " 

-  VCLIP[(^ 

-  fi)  mod  2e,  -co,  co]  +  $ 

(6-4) 

■  ■  . 

'b  ;  a  <  b 

VCLIP (a , 

b,  c)  «  - 

a  ,  b  i  a  i  c 

(6  5) 

,c  ,  a  >  c 

is  a.  Numerix  array  processor  math  library  firtiCtlbri.  (A  Convenient  way 
to  perform  the  modulo  2r  function  Is1  by' 'the'  successive  operations  RECO 
and  POLAR,  which  converts  the  moduluS/itf®5©  tb  real'/imaglnary  and  then 
back  to  modulus/phase*)  This  Fourier  dbniain  bpeidtiort  Is  illustrated 
in  Figure  6-2. 
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Phase 


y  +  <(>t  =  Degraded  phase 

<t>  =  Computed  phase 
8  =  co  -  Phase  limit 

or  =  Phase  error  standard  deviation 


FIGURE  6-2.  PHASE  VARIANCE  ALGORITHM.  Th«  plUM  t'M  1*  C0Mtr»1n«d 
to  lie  within  ca  of  the  given  phase,  #(u). 
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Several  variations  of  the  phase  variance  algorithm  were  attempted 
for  both  the  cases  of  a  real,  nonnegative  object  and  a  complex-valued 
object.  In  the  object  domain  we  can  use  either  the  hybrid-output  (HIO) 
or  error-reduction  (ER)  algorithm  while  one  employs  the  phase  variance 
algorithm  in  the  Fourier' domain.  We  refer  to  these  two  combinations  as 
PVHlb  and  PVE'r, *  respectively,.  Although  HIO  usually  outperforms  ER,  we, < 

found  that  PVER  usually  outperforms  PVHIO. 

.  4  . 

’  i 

Several  values  of  the  variance  factor  c  were  tested.  The  value  of 

•  ‘  i  ' 

c  should  be  small  enough  to  reinforce  the  given  phase  values,  but  large 
enough  to  allow  the  phase  the  freedom  to  adjust  to  become  more 
consistent  with  the  more  accurate  Fourier  modulus  data.  Generally  c  In 
the  range  of  0.6  to  1.0  worked  the  best.  Increasing  or  decreasing  the 
value  of  c  as  the  Iterations  progressed  did  not  seem  to  Improve 
convergence. 

Two  different  Initial  estimates  were  tested.  One  was  the  Image 
gQ(x)  obtained  using  the  noisy  phase  estimate,.  This  is  equivalent  ,tp 
starting  In  the  Fourier  domain  with  phase  f  •  $,  the  noisy  phase.  The 
second  was  an  Image  consisting  of  the  support  constraint  filled  with 
uniformly  distributed  random  numbers.  Host  often  the  nolsy-phase 
Initial  estimate  performed  better  than  the  random  Initial  estimate. 

i  *  .  ■  *■-••• 

It  was  found  that  the  phase  variance  algorithm  would  Improve  the 
estimate  for  several  Iterations,  and  then  it  would  stagnate.  The 
reason  for  the  stagnation  appears  to  be  that  the  outlying  noisy  phase 
values  for  which  the  phase  error  l^el  >  c o  are  inconsistent  with  the 
phase  variance  constraint  [Eq.  (6-3)].  We  found  It  best  to  sUp 
enforcing  the  phase  variance  constraint  at  this  point  and  thereafter  to 
only  enforce  the  Fourier  modulus  constraint  (along  with  the  Image- 
domain  support  constraint).  That  Is,  after  the  phase  variance 
algorithm  stagnates,  continue  with  the  traditional  iterative  transform 
algorithm,  using  cycles  of  HIO  and  ER. 
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After  testing  numerous  combinations  of  algorithm  types  and 
algorithm  parameters,  we  arrived  at  the  following  combination  that 
seemed  to  work  the  best  on  average.  Perform  twenty  Iterations  of  PVER 
with  c  ■  0.8,  then  ten  Iterations  of  ER,  and  finally  several  cycles, 
each  cycle  consisting  of  20  HIO  (beta  ■  0.7)  and  10  ER,  until 
stagnation  (no  further  progress)  occurs.  After  every  other  cycle, 
enlarge  the  support  constraint  by  adding  to  It  each  nearest-neighbor 
pixel  that  was  previously  outside  the  support.  In  order  to  reduce 
sldelobes  In  the  Image,  to  make  the  support  constraint  more  effective 
when  diffraction  effects  are  Included,  the  Fourier  modulus  should  be 
weighted  with  an  apodlzlng  function.  For  the  experiments  described  in 
what  follows,  we  used  a  weighting  function  proportional  to  the 
autocorrelation  of  a  circle  [giving  an  Impulse  response  of  the  form 
(Jj(r)/r)  to  the  complex  Image]. 

The  progress  of  the  Iterative  transform  algorithm  Is  monitored  by 
computing  the  object  domain  error  metric, 


ODEM 


Z  lo'Mi* 

m. 


Z  lf(x)lz 


(6-6) 


where  7  is  the  set  of  points  at  which  g'(x)  violates  the  object-domain 
constraints.  It  Is  a  measure  of  how  close  the  output  Image,  g'(x),  Is 
to  satisfying  the  object-domain  constraints.  For  these  digital 
simulation  experiments,  In  which  we  also  know  the  actual  object,  we  can 
also  compute  the  absolute  error, 


ABSERR  - 


Z  ig'fx  -  x.)  -  f(x)iz 


Z  if(*) iz 


(6-7) 


119 


where  x_  Is  the  shift  of  the  output  Image  g 4 (x)  that  maximizes  jts 

3  •  >  •'  *  *  r 

.correlation  with  the  true  object  f(x).  For  Images  that  are 
recognizable  and  have  some  utility,  ABSERR  Is  typically  below  0.5.  For 
images  that  are  good  representations  of  the  object,  ABSERR  Is  typically 
about  0.3  or  less. 

Figure  6-3  show  ABSERR  as  a  function  of  Iteration  number  for  the 
case  of  using  the  original  iterative  transform  algorithm  with  a.  random 
initial  estimate,  for  the  case  of  a  real-valued;  nonnegative  object 
when  there  Is  no  Fourier  phase  Information.  Each  curve  represents  a 
different  trial  of  the  algorithm  with  a  different  random'  start.  F6r 
the  majority  of  the  cases  the  HIO  algorithm  converges  to  a  good  Image, 
whereas  the  ER  algorithm  .  rarely  does  for  real,  nonnegative  objects. 
When  the  algorithm  was  started  with  a  nblsy  phase  estimate,  without 
reinforcing  It  during  the  Iterations,  it  did  not  Improve  the 
performance  significantly. 

Phase  errors  used  for  these  experiments  were  generated  using 
McGlammery's  method  [6.3].  These  phase  errors  are  similar  to  those 
that  would  result  from  atmospheric  turbulence.  The  adjustable 
parameters  of  the  phase  errors  are  the  standard  deviation,  <r,  and  the 
correlation  length,  corl. 

.Figure  6-4  shows  the  convergence  (In  terms  of  ABSERR)  of  the  ph.a^e 
variance  algorithm  for  several  values  of  the  variance  factor,  c.  From 

i  -  ">  • 

this  we  see  that  the  optimum  value  of  c  for  this  case  is  about  0.8  to 
1.0.  Figure  6-5  shows  the  same  thing  In  terms  of  ODEMl  Although  the 
values  of  ODEM  are  typically  much  less  than  the  values  of  ABSERR,  they 
correlate  fairly  well  with  the  values  of  ABSERR. 

Figure  6-6  show  blurred  Images  and  images  reconstructed  by  the 
phase  variance  algorithm  for  three  different  phase  errors  for  the  case 
of  a  real,  nonnegative  object.  In  all  three  cases,  with  no  noise  on 
the  Fourier  modulus,  the  quality  of  the  reconstructed  images  was  the 
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Error  in  Reconstruction  (ABSERR) 


Reconstruction  Using  the  Original  Iterative  Algorithm 

Random  Inputs 


FIGURE  6-3.  CONVERGENCE  OF  THE  STANDARD  ALGORITHMS. 
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Error  in  Reconstruction  (ABSERR) 


Reconstruction  Using  Phase  Variance  Algorithm 


&  =  2tt/4,  corl=6 


FIGURE  6-4.  CONVERGENCE  (ABSERR)  OF  THE  PHASE  VARIANCE  ALGORITHM  FOR 
VARIOUS  C. 
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Error  in  Reconstruction  (ODEM) 


Reconstruction  Using  Phase  Variance  Algorithm 

<r  —  2rr/4,  corl^ 


FIGURE  6-5.  CONVERGENCE  (ODEM)  OF  THE  PHASE  VARIANCE  ALGORITHM  FOR 
VARIOUS  C. 
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FIGURE  6-6.  RECONSTRUCTION  OF  REAL,  NONNEGATIVE  IMAGES  BY  THE  PHASE 
VARIANCE  ALGORITHM,  (a)  Object;  (b)  support  constrain;  (c) -(d)  blurred 
images,  with  (c)  phase  errors  a  •  */2  radians  and  corl  ■  6  pixels,  with 
d)  a  “  *72  and  corl  «  30,  and  with  (e)  a  -  *75  and  corl  «  6;  and  (f)  - 
(h)  corresponding  images  reconstructed  by  the  phase  variance  algorithm. 
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same  in  all  three  cases,  although  the  convergence  was  faster  for  the 
cases  with  the  smaller  given  phase  error. 

The  reconstruction  of  complex-valued  objects  Is  much  more  difficult 
since  one  no  longer  has  the  powerful  nonnegativity  constraint.  For  the 
complex-valued,  speckled  versions  of  the  object  shown  In  Figure  6-6, 
the  conventional  Iterative  transform  algorithm,  when  starting  with  a 
random  Initial  estimate,  failed  to  reconstruct  a  recognizable  Image. 
When  the  noisy  phase  was  used  to  start  the  algorithm,  however,  the 
conventional  Iterative  transform  algorithm  Improved  the  Image  quality 
substantially,  although  the  reconstructed  Image  remained  Imperfect. 
The  phase  variance  algorithm  similarly  reconstructed  a  substantially 
Improved,  but  Imperfect  Image.  The  results  from  the  phase  variance 
algorithm  were  slightly  better  than  those  of  the  conventional  algorithm 
for  these  cases.  This  Is  Illustrated  by  Figures  6-7  to  6-9. 

Figure  6-10  shows  the  convergence  for  the  case  of  Fourier  modulus 
data  corrupted  by  photon  noise.  The  Iterations  Improved  the  RMS  error 
of  the  Image  for  the  cases  of  more  than  120  photons  per  aperture-plane 
speckle  (or  10®  total  photons).  However,  for  lower  light  levels,  the 
algorithm  can  make  the  image  worse.  This  happens  when  the  Fourier 
magnitude  data  are  noisier  than  the  phase  data;  then  adjusting  the 
phase  to  be  more  consistent  with  the  noisy  modulus  data  Is 
counterproductive.  In  such  a  case  It  would  actually  make  sense  to 
adjust  the  modulus  to  be  more  consistent  with  the  phase  data.  Figure 
6-11  shows  the  RMS  error  of  the  reconstructed  Image  as  a  function  of 
the  total  number  of  photons.  Figure  6-12  shows  three  of  the 
reconstructed  Images.  3  x  10®  total  photons  (120, photons  per  speckle) 
were  required  to  obtain  Improved  Imagery. 

In  summary,  we  have  developed  a  new  variation,  the  phase  variance 
algorithm,  of  the  Iterative  transform  algorithm  which  reconstructs  a 
fine  resolution  Image  when  degraded  Fourier  phase  data  is  available. 
For  real,  nonnegative  objects  It  converges  faster  and  more  reliably 
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Reconstruction  using  Iterative  Phase  Retrieval  Algorithm 
No  Phase  Information  vs.  Noisy  Phase  Input 


FIGURE  6-7.  CONVERGENCE  OF  THE  CONVENTIONAL  ALGORITHM  FOR  A  COMPLEX¬ 
VALUED  OBJECT. 
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Error  in  Reconstruction  (ABSERR) 


Reconstruction  using  Phase  Variance  Algorithm 
No  Phase  Information  (Reg.  recon.)  vs. 
Noisy  Phase  Input,  Phase  Variance  Algorithm 


Number  of  Iterations 


FIGURE  6-8.  CONVERGENCE  OF  THE  PHASE  VARIANCE  ALGORITHM  FOR  A  COMPLEX- 
VALUED  OBJECT. 
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RECO^TPUCtON  OF  LOMPLEX-iiftLiJED  OB.IECt  USING  PHASE  UWIf^tCE  AlGOPITHM 


u,'  PF.irriN.  image  'H'  image  -i>  i»ec  on.  ihhge 


FIGURE  6-9.  RECONSTRUCTION  OF  COMPLEX-VALUED  IMAGES  BY  THE  PHASE 
VARIANCE  ALGORITHM.  (a)  Object;  (b)  support  constraint;  (c)  image 
reconstructed  with  no  phase  information;  (d) -(f)  blurred  imaaes,  with 
(d)  phase  errors  a  -  ir/2  radians  and  corl  =  30  pixels,  with  (e)  a  u  rt 5 
and  corl  *  6,  and  with  (f)  a  =  *75  and  corl  »  30;  and  (g)  -  (i) 
corresponding  images  reconstructed  by  the  phase  variance  algorithm. 
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Reconstruction  Using  Phase  Variance  Algorithm 

<r  =  27t/10,  corl  =  6 

Poisson  noise  added  to  Fourier  magnitude 


Number  of  Iterations 


FIGURE  6-10.  CONVERGENCE  OF  THE  PHASE  VARIANCE  ALGORITHM  FOR  A 
COMPLEX-VALUED  IMAGE  FROM  NOISY  FOURIER  MODULUS  DATA.  For  a  *  r/5  and 
corl  •  6. 
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Reconstruction  using  Phase  Variance  Algorithm 
Poisson  noise  added  to  Fourier  intensity 
a  —  2tt/T0,  corl  =  6 


FIGURE  6-11.  RMS  ERROR  OF  THE  RECONSTRUCTED  IMAGES  AS  A  FUNCTION  OF 
THE  TOTAL  NUMBER  OF  PHOTONS. 
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FIGURE  6-12.  RECONSTRUCTION  OF  COMPLEX-VALUED  IMAGES  BY  THE  PHASE 
VARIANCE  ALGORITHM,  (a)  -  (c)  Noisy  Fourier  modulus  estimates  (their 
squares,  the  intensities,  were  subjected  to  photon  noise),  with  (a) 

107,  (b)  106,  and  (c)  3  x  IQ5  total  photons;  (d)  -  (f)  images  degraded 
by  the  phase  error;  and  (g)  -  (i)  the  corresponding  reconstructed 
images. 
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than  the  conventional  Iterative  transform  algorithm.  For  complex¬ 
valued  objects,  which  are  difficult  to  reconstruct,  it  reconstructs 
images  of  quality  substantially  better  than  that  of  the  blurred  images 
given  by  the  available  noisy  phase  data. 


Other  variations  of  the  phase  variance  algorithm  are  possible  which 


may  yield  improved  performance.  Rather 
new  phase  estimate  such  as  Eq.  (6-3) 
threshold  value,  It  may  be  better  to 
continuously  and  smoothly  with  the  data, 
would  be 


than  using  a  formula  for  the 
which  abruptly  changes  at  a 
have  a  formula  that  changes 
An  example  of  such  a  formula 


-  6  +  ca  ln[l  +  II  -  0!/(c«r)]  s1gn(^  -  0)  (6-8) 

which  is  approximately  equal  to  l  for  1^  -  01  «  c a  and  departs  slowly 
from  the  neighborhood  of  0  when  1^-01  »  c a. 

Another  interesting  possibility  Is  to  use  the  same  type  of 
operation  on  the  modulus  of  the  Fourier  transform.  That  Is,  rather 
than  substituting  the  measured  Fourier  modulus  for  the  computed  Fourier 
modulus,  allow  the  Fourier  modulus  wander  from  the  measured  value 
according  to  the  amount  of  noise  present  in  the  Fourier  modulus  data. 
Such  an  algorithm  would  reconstruct  the  phase  from  the  modulus  or  the 
modulus  from  the  phase  depending  on  which  has  the  higher  slgnal-to- 
noise  ratio  at  any  given  point  in  the  Fourier  domain. 
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7.0  2-D  SHEAR  AVERAGING 


Shear  averaging  is  an  algorithm  invented  at  ERIM  [7.1]  originally 
for  correcting  one-dimensional  (1-D)  phase  errors,  as  occur  in 
synthetic  aperture  radar  (SAR).  Here  we  generalize  It  and  apply  it  to 
the  case  of  2-D  phase  errors  as  would  be  encountered  in  imaging 
satellites  as  described  in  Section  6.  In  what  follows  it  is  seen  that 
a  2-D  extension  to  shear  averaging  is  feasible  if  the  phase  errors  are 
slowly  varying. 

7.1  2-D  SHEAR  AVERAGING  THEORY 

7.1.1  The  2-D  Shear  Averaging  Algorithm 

As  in  previous  work  [7.1],  let  the  ideal  Fourier  transform  be 

F(u,  v)  »  IF(u,  v) I  exp[il(u,  v)]  (7-1) 

and  the  actual,  measured  Fourier  transform,  with  phase  errors  f)#(u,  v) 
be 


G(u,  v)  •  F(u,  v)  exp[i#e(u,  v)]  (7-2) 

where  u  ■  0,  1,  ....  NQ  -  1  and  v  ■  0,  1,  ...,  MQ  -  1.  For  these  2-D 
phase  errors  we  form  the  shear  average 

Soa(u,  v)  -  Z  G(u\  V)  G*(u',  V  -  a)  (7-3) 

oa  (u'.vTeB^ 

where  Buv  is  a  set  of  points  (u\  v‘)  (the  region  of  summation) 
centered  at  (u,  v),  and  a  is  a  lag  smaller  than  the  speckle  size 
(correlation  length)  of  F.  More  generally  a  weighted  summation  can  be 
performed.  Henceforth  the  symbol  B  under  the  summation  means  (u 1 ,  v') 
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In  the  earlier  work  [7.1]  it  was  assumed  that  the  phase  error  is 
one-dimens tonal  (1-D),  i.e.,  ^e(u,v)  ■  ^0(v).  Then  one  can  sum  over 
one  entire  line,  i.e.,  Buv  consists  of  (u',v),  u'  -  1,  2,  ...,  NQ, 
where  NQ  is  the  number  of  samples  in  the  u-dimension.  Then  SQa(u,v)  - 
Sa(v)  is  a  function  of  v  only  and  its  phase  can  be  summed  to  estimate 

a 

4  (v)  [7.1].  The  next  most  complicated  case  is  for  the  phase  error  to 
be  2-D  but  separable,  i.e.,  ^e(u,v)  ■  #eu(u)  +  i*ev(v).  If  one  performs 
the  summation  of  Eq.  (7-1)  again  over  one  entire  line,  one  again  gets 
the  same  result  as  for  the  1-D  cases  $„_(u,v)  ■  S,(v)  is  a  function  of 
v  only  and  its  phase  can  be  summed  to  estimate  fev(v).  Similarly  the 
phase  of  Sj30(u)  [see  Eq.  (7-8)]  can  be  summed  to  estimate  ^eii (1‘)  • 
Consequently  the  separable  case  can  easily  by  handled  as  two  1-D 
problem  with  errors  in  each  dimension  the  same  as  for  the  1-D  case 
[7.1].  In  what  follows  we  consider  the  fully  2-D  case. 

If  the  fully  2-D  phase  error  ^e(u,v)  is  smoothly  varying,  then  we 
can  consider  a  region  of  summation  Buy  having  an  area  over  which  can 
be  approximated  by  a  Taylor-series  expansion  Including  only  linear 
terms: 


Mu''  v')  "  ^u*  v)  +  c10(u’  v)  (u‘  "  u>  +  c01*u'  *v'  "  ^  *  *7"4* 

(Later  we  will  consider  the  effect  of  higher-order  terms.)  Inserting 
Eqs.  (7-2)  and  (7-4)  into  Eq.  (7-3)  and  simplifying  yields 


soa(u,  v)  ■  exp  [la  cQ1(u,  v)]  ^  F(u',  V)  F*(u\  V  -  a)  .  (7-5) 


As  in  the  1-D  case  [7.1],  provided  that  the  area  of  summation  is  large 
enough,  we  can  approximate  the  summation  by  the  ensemble  average 

J(0,  a)  ■  1^(0,  a)  *  <F(u' ,  V)  F*(u\  v'  -  a)>  (7-6) 


138 


where  J(Au,  Av) ,  the  mutual  intensity,  Is  the  Fourier  transform  of  the 
underlying  incoherent  image.  Then  we  have 

SQa(u,  v)  ~  exp[ia  cQ1(u,  v)]  I  /i(0,  a)  .  (7-7) 

Similarly 

SjjQ(u,  v)  ■  ^  G(u 1 ,  V)  G*(u'  -  b,  V) 

“  exp[ib  C|q(u,  v)]  I  /i(b,  0)  .  (7-8) 

Then  the  phases  of  SQa  and  Sbo  are 

eoa(u-  v)  “  «  cQ1(u,  v)  +  arg[^(0,  a))  +  2jr  poa(u,  v)  (7-9) 
and 

0bo(u,  v)  "  b  cio^u»  v)  +  ar9U(b,  0)]  +  Zw  pb0(u,  v)  (7-10) 

respectively,  where  pga  and  pbg  are  integers  that  allow  for  the  fact 
that  the  phase  is  computed  modulo  Zr.  If  the  values  of  a  and  b  are 
chosen  to  be  small  compared  with  the  correlation  distances  of  both  p 
and  j  ,  then  all  the  terms  in  Eqs.  (7-9)  and  (7-10)  will  be  small  and 
Poa  *  pbo  ■  0.  In  the  analysis  that  follows  we  will  not  make  this 
assumption. 

Since  $oa(u,  v)  and  ®b0(u,  v)  represent  phase  derivatives,  by 
integration  over  0oa(u,  v)  and  0bo(u,  v)  one  can  arrive  at  an  estimate 
of  f*e(u,  v).  This  could  be  done,  for  example,  by  first  integrating  in 
the  v  direction  for  a  fixed  u,  then  integrating  in  the  u  direction  for 
each  value  of  v.  The  geometry  and  spacings  of  the  regions  of  summation 
Buy  can  take  several  forms.  In  what  Immediately  follows  we,  give  a 
simple  generic  form  that  lacks  detailed  accuracy  but  explains  the 
principal.  For  example,  suppose  Buy  is  a  rectangular  area  centered  at 
(u,  v)  of  length  mg  in  the  v  direction  and  n0  in  the  u  direction.  Then 
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we  could  compute  on  a  grid  with  spaclngs  mB  and  n0  by  first  summing 
in  the  v  direction  along  u  ■  0: 


(0,  0)  s  0  (7-11*) 

®ia(°'  "V  “  *oa(0'  m 

*  %  £  <^(0,  m'mB)  +  m(mB/a)  arg[/*(0#  a)) 

m 

+  (mB/a)  2t  IJ  poa(0,  *  (Mib) 

Next,  for  each  v  ■  mmB  we  sum  over  u: 

?e(nnB,  mraB)  -  0^(0,  mmB)  +  (nB/b)  ^  *bo(n'nB,  mmB)  .  (7-12a) 

m  n 

-  *B  c01<0'  "'"B*  +  "B  £,  C10<"‘V  "V 

+  «i(mB/ a)  arg[/t(0,  «)]  +  n(nB/a)  «rg|ji(b,  0)] 

m 

+  (mB/a)  2sr  T.  Poa(0,  m'n»B) 
nr-1 

n 

+  (nB/b)  Zw  jjST  Pbo(n'nB,  mmB)  .  (7-12b) 

As  long  as  (mg/a)  and  (ng/b)  are  integers,  then  the  last  two-  summations 
add  integer  2i r  phase  which  Is  unimportant  and  can  be  Ignored.  The 
m(mB/a)  argfy(0,  a)]  +  n(nB/a)  arg[/j(b,  ,0)]  are  linear  phase  terms 
which  shift  the  Image  but  do  not  cause  blurring,  and  so  they  can  be 
Ignored,  The  first  two  terms  mBEcol(0,  m'mB)  +  nBEc1Q(n'nB,  mmB)  are 
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sums  over  the  derivative  of  the  phase  error  which  should  give  a  good 
approximation  to  as  long  as  Eq.  (7-4)  is  accurate  over  each  Buv. 

/\ 

Two  sources  of  error  cause  inaccuracies  in  (1)  statistical 
errors  in  approximating  the  finite  sum  in  Eqs.  (7-7)  and  (7-8)  by  I p 
and  (2)  phase  errors  ^  that  have  higher-order  terms  within  the  region 

V 

of  summation. 

7.1.2  Residual  Phase  Error  Due  to  Statistical  Error 

The  residual  phase  error  due  to  the  approximation  of  the  summation 
over  the  product  of  the  F's  by  the  ensemble  average,  which  is  the  error 
of  arg[/i(0,  a)]  and  arg[/i(b,  0)],  is  similar  to  that  in  the  1-D  case 
[7.1]  and  has  standard  deviation  given  by 

a - I - 

a  JZRjp  l/i(0,  a) I 

and 

_ _ 1 . 

b  iptb,  o)i 

where  in  this  case  Ng  is  the  number  of  independent  samples  of  F  in  the 
region  B.  For  region  B  of  ng  by  mg  samples,  Ng  ■  (ng/nc)  (mg/mc) , 
where  n  by  m  is  the  size  of  an  independent  sample  of  F(u,v).  The 

w  V* 

estimate  of  the  phase  error  across  the  width  ntg  would  be  (mg/a)  0Q  . 
Hence  the  variance  of  the  phase  error  estimate  across  the  width  mg  is 

(mg/a)z  aZ  «  (mg/az)  ncmc/[(2ngmg) l/i(0,a) I2] 

-  mBncn)c/[2a2nBlA»(°,a) 1 2)  .  (7-14) 


(7-13a) 

(7-13b) 
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If  the  total  array  size  Is  NQ  ■  NnB  by  MQ  *  MmB,  then  the  variance  of 
the  phase  error  at  the  far  edge  of  the  array,  assuming  the  simple 
summation  of  Eq.  (7-11),  would  be 


fM. 


mBj  2a^nB  l/i(Q,a)l^  2a^n0  l/*(0,a)l^ 


man_m_ 
o  c  c 


1  “ 


M.n.rtL 
_ o  c  c 


(7-15) 


where  It  Is  assumed  that  the  H  ■  (M0/mB)  errors  are  uncorrelated  over 
the  sum. 


A  similar  result,  exchanging  m's  and  n's,  holds  for  summation  In 
the  orthogonal  dimension,  and  so  the  variance  of  the  error  In  the 
corner  farthest  from  the  beginning  corner  Is 


M.n„nL 
o  c  c 


N_n_m_ 

+  o  c  c 

I"2  ni.2  2  i .. /k  1 2 


2a*ng  l/i(0,a)r  2b*  m0  l/i(b,0) 


for  the  simple  summation  approach. 


(7-16) 


Consider  the  case  of  nc  ■  mc  “  2  samples  and  a  ■  1.  Then,  from  Eq. 
(7-15)  we  see  that  unless  n0  Is  comparable  to  MQ|  will  be 
unacceptably  large  (much  greater  than  one  radian). 

A  A 

However,  there  are  multiple  paths  to  sum  from  ^e(0,  0)  to  #e(NnB, 
MmB).  Methods  used  to  reconstruct  2-D  phase  functions  from  phase 
differences  (least-squares  solutions,  for  example)  can  be  used  here; 
then  the  variance  of  the  residual  phase  error  should  be  much  less  than 
that  given  by  Eq.  (7-16).  Since  2-D  least-squares  methods  reportedly 
yield  a  phase-error  estimate  that  has  an  error  comparable  to  the  error 
In  a  single  phase-error  difference  estimate,  the  variance  of  the  2-D 
phase  error  estimate  should  be  similar  to  that  of  Eq.  (7-14). 
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Note  from  Eq.  (7-15)  that  the  statistical  error  of  the  simple  1-D 
summation  is  independent  of  mg.  Therefore,  as  far  as  that  error  is 
concerned,  the  width  of  Buy  can  be  anything  one  desires.  Since  using  a 
narrow  B  will  reduce  the  effects  of  nonlinearities  in  0  (which  are 
analyzed  later),  it  seems  that  to  compute  $oa(u,v)  one  would  want  to 
use  Buy  of  width  one  sample  in  the  v  dimension  by  nB  samples  in  the  u 
dimension,  where  n0  is  the  greatest  length  that  does  not  run  into 
severe  nonlinearity  problems.  On  the  other  hand,  to  compute  S|30 (u , v) 
one  would  want  Buv  to  be  of  width  mg  samples  in  the  v  dimension  (the 
largest  nig  that  avoids  severe  nonlinearity  problems)  by  one  sample  in 
the  u  dimension.  Thus  we  are  lead  to  using  very  different  sets  of 
points  Buv  for  the  summation  of  SQa  and  S^,  and  different  sampling 
grids  would  result  as  well. 

In  order  to  arrive  at  a  rectangular  grid  from  which  we  could 
proceed  with  a  least  squares  solution,  we  could  first  sum  across  blocks 
of  width  mg  of  the  0Qa  and  down  blocks  of  width  n0  of  ®bo  to  get 
samples  on  a  grid  with  spacings  n0  by  mB.  The  variance  of  the  phase 
error  across  a  block  of  width  mB  is  given  by  Eq.  (7-14).  Take  the 
difference  between  the  phase  values  at  the  beginning  and  end  of  each 
block  to  estimate  the  phase  difference  between  those  two  points.  These 
would  have  the  same  variances  as  mentioned  above.  Then  we  could 
proceed  with  the  least  squares  solution. 

However,  from  Eq.  (7-14)  we  see  that  when  summing  over  just  a 
single  block  of  width  mB,  the  variance  of  the  error  is  proportional  to 
mB<  Therefore  it  would  be  best,  when  using  a  least-squares  summation, 
to  sum  over  narrow  blocks  to  reduce  the  statistical  error  as  well  as 
the  nonlinearity  error.  This  would  suggest  summing  over  blocks  of 
width  mc  for  0Qa  and  of  width  nc  for  0^.  If  the  normal  assumptions 
about  the  least-squares  phase  reconstruction  were  true,  then  the  2-D 
phase  estimate  variances  would  be 
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M2  °\ 


ncm^/[2a2n^  i/*(0,a)l2] 


(7- 17a) 


for  the  v  dimension  for  which  Ng  *  (nB/nc)  x  1*  and 

(nc/bK  '  ncrnc/[2t>2n’B  'MM)!2]  (7-  17b ) 

for  the  u  dimension,  for  which  Ng  -  1  x  (mB/mc).  However,  this  ignores 
the  fact  that  the  error  of  0oa(u,v)  is  correlated  over  ng  samples  in 
the  u  dimension  and  0|jO(u,v)  is  correlated  over  mg  samples  in  the  v 
dimension.  Consequently,  the  variance  of  the  2-D  phase  error  estimate 
may  be  closer  to  ncmc/[2a  l^i(0,a)l  ].  A  further  refinement  of  this 
analysis  will  be  necessary  to  arrive  at  a  more  precise  statement  of  the 
residual  phase  error  for  the  2-D  case. 

7.1.3  Convolutional  Processing 

An  alternative  processing  scheme  is  suggested  by  Eq.  (7-3),  which 
is  essentially  a  convolution  of  G(u',  v ' )  G  (u',  v'  -  a)  with  B0Q. 
Letting  the  functional  representation  of  Buy  be 

B„>'.  V)  ■  {*  f0r  (U''  v'>  e  V) 
uv  10  for  (u\  v')  k  B(u,  v) 

-  B0Q(u'  -  u,  v'  -  v)  ,  (7-18) 

Eq.  (7-3)  can  be  written 

Soa(u,  v)  =  B00(u'  -  u,  v'  -  v)  G (u ' ,  v')  G*(u\  v'  -a) 

U  ,  V 

-  tG(u,v)  G*(u,  »  -  a)]  8  B00fu,v)  (7-19) 
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where  «  denotes  cross-correlation.  More  generally,  BQ0(u ' , v ' )  could  be 
a  nonbinary  function.  A  cross-correlation  can  be  computed  by  two 
FFT's,  a  product,  and  an  inverse  FFT.  For  Buv  covering  just  a  small 
area,  direct  cross-correlaticn  is  more  efficient  than  the  FFT  method. 
By  this  approach,  then,  one  arrives  at  SQa(u,  v)  for  each  sample  of  (u, 
v)  [and  similarly  S^0(u,  v) ] ,  not  just  for  the  coarse  grid  (nnB,  mmB) . 
Integrating  or  performing  a  least-squares  fit  over  this  fine  sampling 
of  the  phase  derivative  is  an  alternative  to  the  use  of  the  coarser 
grid. 

7.1.4  Residual  Phase  Error  Due  to  Higher-Order  Phase  Errors 

The  residual  phase  errors  given  by  Eq.  (7-14)  consider  only  the 
result  of  averaging  over  a  finite  number  of  pixels  to  estimate  the 
ensemble  average.  A  second  source  of  error  is  the  fact  that  the  phase 
errors  are  not  constant  over  the  area  of  integration.  In  what  follows 
is  analysis  of  that  component  of  the  residual  phase  error. 


Now  consider  phase  errors  of  the  more  general  form 


f*(u',  V)  =  ^  Cjk(u,  v)  (u1  -  u)j  (v’  -  v)k  (7-20) 


where 


C00(U|  v)  a  *e(u’  v) 


(7-21) 


for  the  region  B(u,  v) .  Inserting  Eq.  (7-20)  into  Eq.  (7-3)  yields 

$oa(u,  v)  *  X  G(u ' ,  V)  G*(u' ,  v'  -  a) 

B 
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21  F(u* ,  v')  F  (u‘,  v 1  -  a) 

B 

expi|jT  Cjk(u,  v)  (u'  -  u)^  (v1  -  v)k 

-  ^  Cjk(u#  v)  (u*  -  u)J  (V  -  v  *  a)k 
^  F (u ' ,  v')  F*(u' ,  v'  -  a)  expij^r  Cjk(u,  v)  (u' 

[(V  -  v)k  -  (v*  -  v  -  a)k]j 
]T  F (u 1 ,  v')  F*(u\  v'  -  a)  expi||^  cjk^u*  ^u' 
•  [k(v'  -  v)k-1  a  -  [g]  (v*  -  v)k"2  a2  +  ... 


-  u)' 


-  u)J 


+  (-D1"1  *k]) 

where 

fkl  kl 

InJ  "  nTTT^nTT  * 


(7-22) 

(7-23) 


That  Is,  the  phase  error  term  c j k ( u ,  v)  (u*  -  u)^  (v'  -  v)k  results  In 
phase  terms  in  Soa(u,  v)  that  are  jth  order  in  u'  and  (k  -  l)th,  (k  - 
2)th,  ...  order  in  v'.  In  particular  all  terms  in  that  are  zero- 
order  (constant)  in  v  (k  »  0)  disappear,  and  are  therefore 
inconsequential . 


Specifically,  consider  the  phase  error  terms  through  cubic: 


v')  *  Mu*  +  cio^u'  “ u) +  coi^v’  "  ^  *  c2o^u'  " u^2 

+  cn(u'  -  u)  (V  -  v)  +  c02(v'  -  v)2 


¥ 
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+  c30(u'  -  u)3  +  c21(u'  -  u)2(v‘  -  v) 

+  c12(u '  -  u)  (v'  -  v)2  +  cQ3(v'  -  v)3  (7-24) 

where  Cjk  ■  Cjk(u,  v).  Then 

Soa(u,  v)  -  ]T  F (u ' ,  v')  F*(u\  V  -  a) 

exp  i{cQ1  a  +  cu(u‘  -  u)  a  +  cQ2  [2a(v'  -  v)  -  a2] 

+  c21(u'  -  u)2  a  +  c12(u*  -  u)  [2a(v '  -  v)  -  a2] 

+  cQ3  [3a(v‘  -  v)2  -  3a2(v'  -  v)  +  a3]} 

■  exp  i(c„,  a  -  c02  a2  +  c„3  a3) 

Z!  F(u\  v')  F*(u‘ ,  v'  -  a) 

B 

exp  a(u*  -  u)  +  2  cQ2  a(v'  -  v) 

2 

+  c21  a(u'  -  u) 

+  c12(u'  -  u)  [2a(v‘  -  v)  -  a2] 

+  cQ3 [3a( v '  -  v)2  -  3a2(v*  -  v)]}  .  (7-25) 

Now  make  one  further  assumption.  Suppose  that  when  the  summation 
is  replaced  by  an  ensemble  average  over  the  realizations  of  F,  we  can 
treat  the  phase  error  terms  as  being  statistically  independent  of  the  F 
so  we  can  replace  Eq.  (7-25)  by 
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Soa(u,  v)  z  I  MO,  a)  exp  1(cQ1  a  -  cQ2  a2  +  cQ3  a3) 

]T  exp  l{cn  a(u*  -  u)  +  2  cQ2  a(v'  -  v) 

+  c21  a(u'  -  u)2  +  Cj2(u'  -  u)  [2a(v‘  -  v)  »  a2] 

+  c03[3a(v'  -  v)2  -  3a2(v*  -  v)]}  .  (7-26) 

Now  further  suppose  that  B00(u',  v ' )  is  symmetric  In  u'  and  in  v' 

and  separable  in  u1  and  v'.  Then,  individually,  terms  In  Eq.  (7-26) 
that  are  odd  in  (u*  -  u)  or  (v '  -  v)  will  have  integrals  (sums)  of 

their  imaginary  part  that  will  be  zero  (making  the  Integral  have  zero 

phase)  and  no  undesired  phase  terms  will  result.  Consequently,  of  the 
phase  error  terms  explicitly  shown  in  Eq.  (7-24),  the  single  terms  that 
cause  undesired  phase  errors  In  the  summation  of  Eq.  (7-26)  are  those 
having  coefficients  c21  and  Cq3<  If  we  assume  small  phase  error 
contributions  due  to  these  terms,  we  can  approximate 

^  exp  i[c21  a(u'  -  u)2  +  cQ3  3a(v*  -  v)2] 

~  [l  +  1c21  a(u‘  -  u)2  +  1c03  3a(v'  -  v)2] 

ng/2  mg/2 

2  nBmB  +  iirig  c2l  a  J  u'2  du 1  +  ing  cQ3  3a  j  v'2  dv* 

-ng/2  -mg/2 

•  nB*B  +  '"b  C21  “V12  *  ,nB  ^03  3*  mB/12 

‘  VbI1  +  ,c21  *nB/12  +  k03  aill/4] 

~  nB«iB  e*pi[c21  «nB/12  +  c.j  *m|/4l  .  (7-27) 
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Therefore  for  f>e  with  terms  up  to  cubic  over  B,  the  phase  of  $oa(u,  v) 
is  given  approximately  by  • 


0Qa(u,  -  ar9 (° i  a)]  +  cQi  a  -  cQ2  a2  +  cQ3  a(aZ  +  mg/4] 

+  c21  ang/12  +  2f  poa(u,  v)  (7-28) 

(where  the  c ^  are  functions  of  u  and  v). 

If  we  ignore  the  fact  that  there  are  higher-order  phase  errors  and 
take  this  phase  to  be  due  to  the  linear  component  only,  then  the  phase 
difference  across  B  will  be  taken  to  be 

(m0/a)  0oa(u,  v)  -  (mg/ a)  arg[^(0,  a)]  +  cQ1  m0  -  cQ2  m0  a 

+  c03  +  "|/4)  +  C21  mB  "b/12 

+  (mg/a)  2 w  poa(u,  v)  .  (7-29) 

At  this  point  consider  what  the  actual  phase  difference  across  the 
center  of  B  is  in  the  v  direction.  For  u1  ■  u,  we  have,  from  Eq. 
(7-24), 

v)  +  cQ1(±  mg/2)  +  cQ2(±  mg/2)2  +  cQ3(*  nig/2)3. 

(7-30) 

-  ^e(u,  v  -  mg/2)  ■  cQ1  m0  +  cQ3  mg/8  .  (7-31) 

The  residual  error,  the  difference  between  Eqs.  (29)  and  (31),  ignoring 
(mg/a)  arg[^(0,  a)]  and  the  2rp  terms,  is 


^e(0,  v  t  mg/2)  -  *e(u, 

Therefore 

>e(u,  v  +  mg/2) 
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res.  error  “  -  c0g  mB  a  +  Cqj  mB(a^  +  nig/8)  +  Cg^  mB  n^/12  .  (7-32) 

Of  these  terms,  -CggHga  will  probably  be  the  worst  since  its 
coefficient  will  ordinarily  be  the  largest. 

The  last  source  of  error  that  we  will  consider  here  are  the  odd- 
function  errors  in  Eq.  (7-26)  that  were  dropped  because  they  added 
nothing  to  the  phase-error  estimate.  Their  deleterious  effect  is  to 
reduce  the  magnitude  of  the  summation  over  B,  thereby  reducing  the 
signal -to-noise  ratio.  For  example,  the  first  term  taken  by  itself 
would  yield 


nB/2 

51  exp  [iejjaCu'-u)]  ■  mB  J  exp(1c11au')  du' 

B  -nB/2 

■  t  mBnB  sinc[c1;lanB/(2r)]  (7-33) 

which  would  go  to  zero  for  Cj^anB  ■  2r,  Therefore,  the  contribution  of 
such  terms  to  over  B  must  be  considerably  less  than  2ir  In  order  to 

v 

avoid  a  significant  loss  In  signal-to-noise  ratio. 

To  minimize  the  residual  phase  errors  in  Eq.  (7-32)  due  to  higher- 
order  phase  errors,  we  would  want  to  choose  small  values  for  and  nB, 
whereas  to  minimize  the  residual  phase  errors  due  to  the  statistics 
[see  Eq.  (7-16)]  we  would  want  to  maximize  m0  and  nB.  The  optimum 
trade-off,  which  depends  on  the  spatial  statistics  of  le,  should  be 
determined.  Another  possibility  Is  to  perform  the  phase  error 
correction  recursively.  Depending  upon  how  is  corrected  in-between 
the  coarse  sampling  grid  (i.e.  what  form  of  interpolation  is  used), 
some  of  the  higher-order  terms  may  be  reduced.  Then  a  second  pass  of 
the  algorithm  may  improve  the  result. 
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It  is  also  possible  to  reduce  these  errors  substantially  with  more 
complex  processing.  If  we  use  the  convolutional  method  for  computing 
(u ,  v)  at  the  finer  sampling  grid,  then  with  closely-spaced  samples 
of  Soa(u,  v)  we  can  estimate  the  quadratic  coefficient  Cq2  and 
compensate  for  it  appropriately.  Extensions  to  this  approach  could  be 
used  to  reduce  other  higher-order  terms  as  well,  including  cross-terms 
such  as  c21-  The  area  of  optimally  using  the  data  Sba(u,  v)  to 

estimate  jie(u  v)  Is  probably  a  very  rich  area  in  which  great 

Improvements  could  be  made.  One  should  investigate  least-squares  and 
bi spectrum-1  Ike  approaches,  for  example. 

In  the  derivations  given  here  we  estimated  edge-to-edge  phase 
errors  across  the  regions  B,  but  then  corrected  them  on  a  center-to- 
center  basis.  Therefore  the  correction  equations  need  to  be  modified 
to  account  for  this  effect. 

7.2  COMPUTER  SIMULATION  AND  RECONSTRUCTION  EXPERIMENTS  ,  . 

Based  on  the  theory  presented  in  Section  7.1,  three  types  of 
spatial -frequency  summations  were  implemented,  as  illustrated  in  Figure 
7-1. 


Given  Soa(u,v)  and  Sho(u,v),  estimates  of  phase  error  differences 
from  the  integrations,  the  method  we  used  for  reconstructing  the  phase 
error  was  the  complex  exponential  phase  reconstruction  algorithm  shown 
in  Figure  7-2,  which  is  taken  from  Reference  7.2.  In  that  figure  Pmn 

A 

is  equivalent  to  exp[i^e(m,n)] ,  Dumn  is  equivalent  to  S^u.v),.  and 
Dwmn  Is  equivalent  to  S„(u,v).  First  a  simple  product  (phase 
summation)  is  performed  along  each  of  the  two  axes,  then  the  Interior 
points  are  built  up  recursively  using  a  summation  over  two  paths.  Next 
several  iterations  are  performed.  ,  In  one  iteration,  each  value  is 
replaced  by  a  summation  of  values  taken  from  the  four  nearest 
neighbors.  The  order  of  the  selection  of  the  values  Is  in  an  outward 
spiral:  first  a  clockwise  spiral,  then  a  counterclockwise  spiral.  In 
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2-0  SHEAR  AVERAGING 


INTEGRATION  GEOMETRIES 


M  M 


Case  1 

Integrate  over  1  x  m8. 

Using  each  pixel , 
Reconstruct  M  x  M  array. 

{"B 

1 

G(u, v)  G*(u,v-a) 

M 

ut 

mB 

G(u,v)  G*(u-b, v) 

For  S01(u,r) 

L 

For  Sb0(u.v) 

mB 

Case  2 

Eh 

Integrate  over  m0  x  mg. 

Using  each  pixel , 

Reconstruct  M  x  M  array, 

G(u,v)  G*(u,v-a) 

G(u,v)  G*(UtP,v) 

Case  3 

Integrate  over  n>8  x  mB  (like  Case  2). 

Using  each  mg  x  nigth  pixel, 

Reconstruct  (M/mJ  x  (M/mg)  array  from  (n»B/a)  SQa  and  (mB/a)  SbQ. 
Interpolate  to  M  x  M  array 


FIGURE  7-1.  INTEGRATION  GEOMETRIES  FOR  2-0  SHEAR  AVERAGING. 
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COMPLEX  EXPONENTIAL  PHASE 
RECONSTRUCTION 

i  ,  ■  '  vV  ,  V 

•  •  .  •  I  I 


First  estimate: 


P11-1 

P12-Du11  P11/|P12| 

P22  -  (Du21  P21  4-  Dv12  P12)/|P22| 
etc. 

i  •  * 

Iteration  Example: 

P23-(DU22  P22  +  DV13  P13 

+  DJ, 23  P24  +  D*  23  P33]/|P23| 
etc. 


FIGURE  7-2.  COMPLEX  EXPONENTIAL  PHASE  RECONSTRUCT, OR. 
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Figure  7-2,  the  division  by  IPmnl  on  the  right-hand  side  Is  to  be 
Interpreted  as:  first  compute  the  right-hand  side  without  IPmnl,  then 
divide  It  by  Its  magnitude  to  arrive  at  a  pure-phase  function  (a  phase 
exponential) . 

t 

Figure  7-3  shows  a  case  for  which  the  object  Is  a  delta  function 

* 

and  the  phase  error  Is  that  shown  In  Figure  7-3(a),  with  or  ■  ir/2  and 
corl-30.  For  display  of  phases,  -r  Is  black,  r  Is  white,  and  the  phase 
Is  wrapped  (modulo  2f).  The  Impulse  response  for  this  phase  error  Is 
shown  In  Figure  7-3(g).  The  phases  and  magnitudes  of  S0fl  and  for 
Case  1  Integration  (see  Figure  7-1)  with  mB  ■  8  are  shown  in  Figures 
7-3(b)-(e).  The  phase-error  estimate,  reconstructed  by  shear  averaging 
In  conjunction  with  the  reconstructor  shown  In  Figure  7-2,  Is  shown  In 
Figure  7-3(f).  From  this  It  can  be  seen  that  the  reconstructed  phase 
error  Is  similar  to  a  smoothed  version  of  the  true  phase  error.  The 
smoothing  Is  due  to  the  value  of  nu,  Figure  7-3 (h)  shows  the  Impulse 
response  due  to  the  residual  phase  error  gotten  by  subtracting  the 
estimated  phase  error  from  the  true  phase  error.  From  this  It  Is  seen 
that  subtracting  the  phase  error  estimate  removes  most  of  the  error, 

Figure  7-4 (a)  shows  the  complex-valued  object  used  for  the 
experiments  that  follow.  Figures  7-4(b),  (c)  show  the  magnitude  and 
phase  of  the  Fourier  transform  of  the  object.  Figure  7-4 (d)  shows  the 
added  phase  error  (cr  ■  r/2  and  corl  ■  30).  Figure  7-4(e)  shows  the 
given  noisy  phase  [(c)  plus  (d),  modulo  2r] ,  and  Figure  7-4(f)  shows 
the  blurred  image  obtained  using  the  noisy  phase. 

Figure  7-5  shows  a  reconstruction  experiment  similar  to  that  shown  ♦ 
in  Figure  7-2,  with  a  case-1  Integration  (see  Figure  7-1)  with  mB  ■  8. 

Comparison  of  the  original  object  In  (a)  with  the  blurred  Image  In  (c) 
and  the  reconstructed  Image  in  (h)  shows  that  2-0  shear  averaging 
corrected  some  of  the  phase  error,  but  left  a  large  residual  phase 
error.  Figure  7-6  shows  the  same  thing  for  mB  ■  32,  and  Figure  7-7 
shows  the  same  thing  for  a  case-2  Integration  (see  Figure  7-1)  for  mB  ■ 
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FIGURE  7-3.  PHASE  ERROR  RECONSTRUCTION 
POINT-SOURCE  OBJECT,  (a)  Phase  error 
phase  of  SbQ,  S  ;  (d) ,  (e)  magnitude 

phase;  (g)  TmpuTse  response  from  (a); 
minus  (f). 


BY  2-D  SHEAR  AVERAGING  FOR  A 
function  (modulo  2*);  (b) ,  (c) 
of  Sl  ,  S  ;  (f)  reconstructed 
(h)  impuTse  response  from  (a) 
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FIGURE  7-4.  DATA  USED  IN  2-D  SHEAR  AVERAGING  RECONSTRUCTION 
EXPERIMENTS.  (a)  Object;  (b) ,  (c)  object's  Fourier  magnitude  and 
phase;  (d)  phase  error;  (e)  noisy  phase  estimate  [(c)  plus  (d) ,  modulo 
2*];  (f)  image  blurred  by  the  phase  error. 
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FIGURE  7-5.  PHASE  ERROR  CORRECTION  BY  2-D  SHEAR  AVERAGING.  Case-1 
integration  and  «  8  were  used,  (a)  The  object;  (b)  the  phase  error; 
(c)  the  blurred  image;  (d) ,  (e)  the  phase  of  Sh  ,  S  ;  (f) ,  (g)  the 
magnitude  of  Sb  ,  S  ;  (h)  the  image  reconstructed°after  2-D  shear 
averaging.  a 
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FIGURE  7-6.  PHASE  ERROR  CORRECTION  BY  2-D  SHEAR  AVERAGING.  Same  as 
Figure  7-5,  except  nig  =  32  was  used. 
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FIGURE  7-7.  PHASE  ERROR  CORRECTION  BY  2-D  SHEAR  AVERAGING.  Same  as 
Figure  7-5,  except  case-2  integration  and  ntg  =  32  was  used. 
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16.  Of  these,  the  case-1  integration  with  mg  *  32  appears  to  yield  the 

best  image  for  this  example.  Larger  values  of  m0  cause  the  integration 

to  be  over  an  area  over  which  the  phase  error  varies  too  wildly, 

whereas  smaller  values  of  nig  cause  the  integration  to  be  over  a  smaller 

region,  increasing  the  statistical  error.  Comparing  Figure  7-  5(c)  ; 

with  7-5(h)  shows  that  2-D  shear  averaging  improves  the  quality  of  the 

image  substantially,  but  far  from  perfectly.  >. 
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8.0  SPACE  OBJECT  IMAGING  SENSORS 


This  section  describes  a  first  rough  cut  at  comparing  the  numerous 
potential  sensors  for  imaging  space  objects  in  earth  orbit  from  the 
ground  using  short  (visible  or  near-IR)  wavelengths.  Because  no 
extensive  investigations  were  performed  to  compare  the  various  imaging 
approaches,  what  is  contained  iri  this  section  should  not  be  considered 
to  be  a  recommendation  of  one  approach  over  another;  rather,  this 
should  be  viewed  as  an  off-the-cuff  listing  of  attributes  and  as  only  a 
first  step  toward  comparing  the  various  approaches.  There  is  a  need  to 
perforin  a  thorough  analysis  comparing  these  numerous  candidate  systems; 
this  was  beyond  the  scope  of  the  present  program,  but  it  is  recommended 
that  such  an  analysis  be  performed  to  establish  the  basis  for 
development  of  future  fine-resolution  imaging  systems.  This  comparison 
is  done  primarily  by  means  of  the  three  matrix  charts  shown  in  Figures 
8-1  to  8-3.  Figure  8-1  covers  the  case  Pf  using  only  laser 
illumination,  Figure  8-2  covers  the  case  of  using  only  incoherent 
illumination  (or  emissive  objects),  and  Figure  8-3  covers  mixed- 
coherence  (coherent/ incoherent)  methods  and  other  miscellaneous 
approaches.  Further  additions  to  the  matrix  could  be  made. 

The  Near-Team  Feasibility  column  indicates  our  opinion  of  the 
feasibility  of  performing  a  successful  experiment  with  present-day 
technology  using  an  existing  single-aperture  telescope  within  the  next 
six  months.  A  successful  experiment  would  be  one  in  which  the 
resolution  of  the  reconstructed  image  is  several  times  better  than  (X 
R/rQ)  at  the  object  without  the  use  of  adaptive  optics,  where  X  ■ 
wavelength,  R  -  range  to  targt  '  and  rQ  ■  Fried's  parameter  (~10cm) . 
The  Large  Distributed  Aperture  column  comments  on  the  difficulty  of 
putting  together  electro-optical  hardware  for  a  large  distributed  array 
of  apertures  suitable  for  imaging  geosynchronous  objects. 
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5.0.1.  SQBOftS  «SUC  CONCKM  USU  liXMIMTUi  W.T 
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FIGURE  8-1.  S.O'I.  SENSORS  USING  COHERENT  LASER  ILLUMINATION  ONLY 


FIGURE  8-2.  S.O.I.  SENSORS  USING  INCOHERENT  ILLUMINATION  ONLY 


S-O.I 


8.1  INCOHERENT-ONLY  SENSORS 

Several  incoherent  approaches,  including  both  aperture-plane  and 
focal-plane  (astronomical  speckle)  interferometry  are  feasible.  Of 
these,  the  one  farthest  along  in  development  is  speckle  interferometry 
(13)  using  either  Knox-Thompson  or  triple  correlation  to  obtain  an 
initial  image  (or  Fourier  phase)  estimate  which  is  refined  by  the 
iterative  transform  algorithm.  This  would  yield  an  Incoherent  image  of 
the  object.  It  requires  hundreds  to  thousands  of  frames  with  different 
realizations  of  atmospheric  turbulence,  and  requires  measurements  on  a 
reference  star  through  an  atmosphere  having  the  same  statistics  as  the 
atmosphere  through  which  the  object  is  Imaged.  It  is  restricted  to 
pre-dawn  or  post-dusk  imaging  while  the  sensor  is  in  night-time  but  the 
object  is  sun  illuminated. 

An  extension  of  this  method  to  large  distributed  apertures  for 
imaging  geosynchronous  satellites  would  be  difficult  due  to  a 
requirement  of  a  common  focal  plane. 

8.2  COHERENT-ONLY  SENSORS 

For  imaging  with  coherent  laser  illumination  only,  all  the 
approaches  are  risky.  The  least  risky,  in  terms  of  image 
reconstruction,  would  be  a  combination  of  Laser  Dual  Plane  (C5)  with 
Imaging  Correlography  (C3) .  Figure  8-4  shows  a  flowchart  of  the  data 
processing  for  this  combined  approach.  Imaging  Correlography  is  the 
collection  of  multiple  aperture-plane  speckle  intensity  patterns,  power 
spectrum  (or  autocorrelation)  averaging,  and  image  reconstruction  by 
the  iterative  transform  algorithm  to  arrive  at  a  moderate-resolution 
(/?m)  incoherent  image.  The  resolution  is  limited  primarily  by 
statistical  averaging  noise  due  to  a  finite  number  of  realizations  of 
the  speckle  patterns.  The  Laser  Dual  Plane  method  is  the  collection  of 
both  an  aperture-plane  and  a  focal-plane  speckle  intensity  pattern,  and 
processing  by  a  Gerchberg-Saxton  type  algorithm  to  arrive  at  the 


163 


FIGURE  8-4.  DATA  PROCESSING  BLOCK  DIAGRAM  FOR  THE  COMBINED  LASER  DUAL 
PLANE  IMAGING  CORRELOGRAPHY  SENSOR. 


aperture-plane  complex  field.  This  gives  the  field,  F  ■  IFIexp(i|>), 
scattered  by  the  object  times  a  phase  factor,  exp(i^),  due  to 
atmospheric  turbulence;  the  result  therefore  has  the  degraded  phase 
f+^a.  Reconstruction  algorithms  described  in  Sections  6  and  7  would  be 
appropriate  for  correcting  the  phase  error,  ^a. 

If  all  the  aperture-plane  snapshots  are  processed  via  the  Imaging 
Correlography  approach  into  a  moderate-resolution  image,  then  it  should 
be  possible  to  use  that  image,  in  conjunction  with  algorithms  for 
setting  upper  limits  ("locator  sets")  on  the  support  of  the  object  from 
the  support  of  its  autocorrelation,  to  determine  a  reasonably  tight 
support  constraint  on  the  object. 

The  support  constraint  from  Imaging  Correlography  plus  the  degraded 

phase  from  the  Dual-Plane  approach  should  make  diffraction-limited 

resolution  (p^)  coherent  image  reconstruction  from  a  single  snapshot  of 

aperture-plane  Intensity  easier  (although  the  question  of  just  how  easy 

it  would  be  has  not  yet  been  investigated).  A  priori  knowledge  of  a 

support  constraint,  which  might  be  known  for  a  "friendly"  object,  would 

also  make  image  reconstruction  easier.  With  the  collection  of  multiple 

frames,  one  has  the  option  of  choosing  which  snapshot  to  process,  and 

the  selection  of  one  for  which  a  strong  glint  is  present  and  favorably 

positioned  would  make  reconstruction  easier  still.  Selection  of  an 

object  that  is  highly  noncunvex  would  also  help.  Noncoherent  averaging 

-1/2 

of  N  reconstructed  images  would  decrease  the  speckle  contrast  to  N  , 
approximating  an  incoherent  image  for  this  and  all  the  other  coherent 
approaches  that  follow. 

This  approach  would  scale  well  for  large  distributed  apertures 
since  the  two  detection  planes  and  the  wavefront  sensing  could  be  done 
independently  for  each  aperture;  however,  the  reconstruction  of  the 
coherent  images  could  suffer  from  the  sparseness  of  the  array. 


8.3  COMBINED  DUAL-PLANE  AND  INCOHERENT  ATMOSPHERE  SENSING 

The  approach  employing  active  illumination  that  is  most  likely  to 
succeed  is  (Ml),  a  combination  of  the  Laser  Dual-Plane  Sensor  (C5)  with 
an  Incoherent  Atmospheric  Sensor.  The  Laser  Dual-Plane  Sensor,  as 
described  above,  yields  the  degraded  phase  y+ya.  A  wavefront  sensor 
such  as  a  shearing  interferometer,  operating  with  incoherent  light  from 
the  object,  yields  the  atmospheric  phase,  Subtraction  of  the 
atmospheric  phase  from  the  degraded  phase  yields  y,  the  phase  due  to 
the  object.  Then  an  image  is  reconstructed  by  inverse  Fourier 
transformation  of  IFIexp(iy). 

This  combined  method  Involves  fairly  complex  hardware?  two 
detector  planes  for  the  laser  wavelength  plus  a  wavefront  sensor  for 
the  Incoherent  light.  Furthermore  it  requires  both  laser  and 
incoherent  (e.g.  sun)  Illumination.  However,  the  phase  retrieval  part, 
finding  the  aperture-plane  field  by  the  Gerchberg-Saxton  algorithm,  is 
low  risk  and  the  wavefront  sensor  Is  already  in  place  as  part  of  the  Cl 
system  at  AMOS. 

This  approach  would  scale  well  for  large  distributed  apertures; 
however  phase  retrieval  will  be  needed  for  inter-aperture  phase  errors. 

8.4  LASER  FOCAL  PLANE  WITH  INCOHERENT  ATMOSPHERE  SENSOR 

This  method  (M2)  is  the  same  as  the  method  above  (Ml)  except  that 
(A)  the  aperture  plane  detector  is  eliminated  and  (B)  the  ordinary 
circular  aperture  is  masked  to  form  an  asymmetric-shaped  aperture. 
Then .instead  of  using  the  Gerchberg-Saxton  algorithm  to  determine  the 
aperture-plane  field  from  the  aperture  and  focal  plane  intensities,  one 
uses  the  iterative  transform  algorithm  to  determine  the  aperture-plane 
field  from  the  focal  plane  intensity  and  the  aperture-shape  support 
constraint.  The  atmospheric  phase  is  subtracted  and  the  Image  Is 
formed  as  in  (Ml)  above. 


8.5  CONVENTIONAL  LASER  IMAGING  APPROACH 


Although  it  is  unlikely  to  be  practical,  we  have  conceived  of  a 
means  whereby  it  would  be  possible  to  reconstruct  a  diffraction-limited 
>  image  from  a  single  focal-plane  intensity  array  in  coherent  light.  It 

is  like  the  Laser  Conventional  approach  (Cl)  but  with  a  long-icoherence 
*  length  laser  with  the  addition  of  an  aperture  of  special  shape.  First, 

as  in  (M2)  above,  one  uses  the  iterative  transform  algorithm  to 
reconstruct  the  aperture-plane  field  from  the  focal  plane  Intensity  and 
the  aperture-shape  support  constraint.  Then  we  use  the  iterative 
transform  algorithm  to  reconstruct  a  coherent  image  from  the  modulus  of 
the  reconstructed  aperture-plane  field  and  a  support  constraint  on  the 
object.  This  method  would  require  a  much  higher  data  slgnal-to-noise 
ratio  than  method  (Ml)  or  (M2)  above,  since  the  modulus  of  the 
aperture-plane  field  must  be  known  more  accurately  for  the  iterative 
reconstruction  of  the  image  than  if  the  phase  is  known  for  simple 
Fourier-transform  reconstruction  of  the  image.  Nevertheless,  this  two- 
stage  reconstruction  approach  is  theoretically  very  Intriguing. 
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